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CHAPTER 1 


INTRODUCTION AND STATEMENT OF PROBLEM 

The localized flow fields induced around buildings 
and other manmade obstructions by the surface winds h&ve 
long been of interest in structural design, but more 
recently, atmospheric motion influenced by such obstacles 
has become important in the design of airports, launch pads 
and other landing facilities for aircraft and spacecraft. 

Of particular interest, in view of the rapidly increasing 
needs in commercial air transportation, is the development 
of helicopter and V/STOL service in large metropolitan areas. 
Because of the possibility of operating these aircraft from 
the tops of buildings, the problem of analyzing building- 
induced flow fields attains a new importance. 

Many difficulties are involved in operating low speed 
aircraft near buildings such as severe crosswinds, induced 
vortex fields, regions of separated flow, and other unsteady 
flow phenomena which make takeoff and landing extremely 
hazardous. There are two possible approaches to minimizing 
the dangers of these problems: the first is related to 

improvements in the design of the aircraft and operating 
procedures during takeoff and landing; and the second 
attacks the source of the flow disturbances by modifying 
airport design in an effort to reduce operational problems 


However, each of these approaches 


created by the wind (1) . ^ 
requires a detailed knowledge of the flow field and turbu- 
lence imposed on buildings and other structures near the 
ground. 

One method for satisfying this need for quantitative 
descriptions of localized flow fields surrounding aircraft 
landing sites, is to make extensive measurements around 
existing facilities, but the problems of locating sensors in 
the most advantageous locations to measure critical flow 
parameters and the vast amount of data required for such a 
description, make this method highly impractical. Further- 
more, such detailed measurements around existing airports 
are of little value in the design of new structures without 
additional theoretical work to establish similarity relation- 
ships between these flow fields. 

In considering elevated V/STOL ports, interest is 
primarily concerned with the flow around buildings; however, 
at present, investigations are limited to qualitative 
descriptions of the effects of isolated structures and 
groups of buildings on atmospheric motions. A number of 
these descriptions are given in References (1), (2), (3) and 

(4) discussing such difficulties as three-dimensional 
effects, separated and wake flows, induced vortices, the 
effects of shear and interaction of turbulence with steady 

■^Numbers in parentheses refer to similarly numbered 
references in the bibliography. 
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flow phenomena, and other complexities which further empha- 
size the necessity for theoretical approaches to the prob- 
lems of atmospheric motion over surface obstructions. 

In view of the above considerations and the importance 
of such flows to aerodynamic design, the need for an analyt- 
ical method for predicting local atmospheric motions 
influenced by buildings and other surrounding topographical 
features is obvious. 

In formulating a model to completely describe this 
complicated flow situation, the complete equations of motion 
for turbulent flow must be considered. However, practical 
methods for carrying out a solution of such equations are 
limited to numerical approaches which are at present in the 
very early developmental stages (5), (6) and (7). Further- 

more, the extremely high financial costs in computing time 
and the restriction of these solutions to specific problems 
limit the use of this method. 

An alternative approach for approximating atmospheric 
motions over surface obstructions is to extend the concepts 
of boundary layer theory, which have been very successful in 
describing aerodynamic flow over surfaces of small curvature 
where viscous regions are thin. In meteor logical work, the 
boundary layer concept has been applied to wind flow over 
changes in terrain roughness (8) , (9) and (10) and wind 

screens (11) and (12), by introducing the idea of the dis- 
turbed boundary layer . 

Consider a fully developed turbulent shear flow 
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approaching a discontinuity in surface condition such as a 
change in surface roughness, as shown in Figure 1 . Because 
the disturbance is convected downstream and diffuses 
vertically, it is assumed that at great distances upstream 
and large heights downstream of this discontinuity, the flow 
is unaffected by the change in surface roughness. The Ekman 
layer may experience an adjustment in height through the 
scaling relationship h ~ u*/f, where u* is influenced by z Q , 
however, this is not expected to affect the flow appreciably 
in the region of interest. This assumption leads to a region 
of influence surrounding the downstream roughness, outside 
of which the flow is characteristic of approach conditions. 
Very near the ground, however, the boundary layer is governed 
by local surface conditions with the outer flow exerting an 
influence only through the boundary conditions between the 
two regions ( 13 ) . 

Although the disturbed boundary layer concept holds a 
great deal of promise, there is little evidence that this 
approach has been applied to flow over buildings and other 
surface obstructions. Therefore, the purpose of this study 
is to investigate the use of this boundary layer concept in 
approximating the localized flow field induced around a 
surface obstruction by the impinging wind. Although the 
method developed within this study can, in general, be 
applied to flow over any two-dimensional body, the solutions 
given will be limited to semi-elliptical cylinders. 
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Figure 1. Concept of the disturbed boundary layer: 6 (x) is the thickness 

of the internal boundary layer; 5' is the displacement of the outer 
flow field required by continuity. 




CHAPTER II 


ANALYSIS. OF ATMOSPHERIC FLOW OVER SURFACE 
OBSTRUCTIONS BY THE TURBULENT 
BOUNDARY LAYER APPROACH 

In applying the concepts of boundary layer theory to 
atmospheric motions over surface obstructions, the following 
two-dimensional model will be assumed. Consider a uniform 
terrain of infinite extent on which a surface obstruction is 
located in the form of a semi-elliptical cylinder as shown 
in Figure 2. The elliptical contour is chosen primarily 
because its geometry simulates low buildings or hills and 
permits parametric variations of the aspect ratio. In 
addition, the inviscid flow field around such an obstruction, 
required for later analysis, is well established in potential 
flow theory. 

The Nature of the Approaching Wind 

Assuming that the homogeneous terrain extends far 
upstream from the elliptical obstruction, the viscous nature 
of the atmospheric motion over the ground produces gradients 
in the mean wind velocity. For conditions of neutral 
stability experimental evidence (29) and (30) confirms that 
the mean wind velocity in the region of the atmosphere 
near the ground is described by a logarithmic expression 
analogous to aerodynamic measurements of turbulent boundary 
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Figure 2. Basic geometry of the flow over a two-dimensional 
elliptical cylinder of major axis 2a and minor axis 2b. 


layers. For smooth surfaces this logarithmic distribution 
is given by the well known "law of the wall" (14) , 




+ C 


( 1 ) 


where < is the universal von Karman constant, C is another 
empirical constant, and u* is the friction velocity defined 
in terms of the surface shear stress as 


u* = /T 0 /p 


This universal velocity distribution was first 
obtained by Prandtl as a consequence of his classical con- 
cept 'of the mixing length for turbulent shear flows. The 
validity of Equation 1 implies a region of constant shear 
near the wall in which the mixing length is assumed propor- 
tional to the normal distance from the surface approaching 
zero at the smooth surface. (A discussion of the mixing 
length theory will be given in a later section of this 
investigation. ) 

In considering the more general case of a rough sur- 
face, the logarithmic velocity distribution given above, must 
be altered to incorporate roughness into the flow near the 
surface. Roughness is a characteristic of most natural 
terrains which is generally described in terms of a vertical 
length scale, z q . This length scale is related to the 
height of the individual roughness elements, but is also 
strongly dependent on the geometry and distribution of the 
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roughness over the surface. Typical values of z Q are given 
in Table I for different types of surface terrain. 

The surface roughness causes an increase in the 
dissipation of momentum near the ground by larger viscous 
losses and more intense turbulent action. However, these 
effects become important only for surfaces on which the 
characteristic length, z Q , exceeds v/u*. In this case, the 
mixing- length does not vanish near the wall but approaches 
the size of the roughness (15) , and is given by the 
following expression: 


i = k(z + Z Q ) 


( 2 ) 


Assuming that the surface roughness is uniformly distributed 
over the ground in a random manner, this relation for the 
mixing length leads to a logarithmic velocity of the 
following form: 



1 

K 


In 


z + z 


(3) 


Therefore, far upstream from the elliptical obstruction, the 
flow field is described by the above equilibrium equation, 
valid for flow over a uniform, rough surface. 


Governing Equations for Atmospheric Flow 

The governing equations for steady, mean incompress- 
ible flow within the atmospheric boundary layer are well 
established in meteorlogical literature (16) and are stated 
in the following relationships: 
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TABLE I 


i 


ROUGHNESS LENGTHS FOR VARIOUS TYPES OF TERRAIN 


Terrain Type 


Surface Roughness 
Length (z q ) in Meters 

Ice and snow 


10~ 7 - 7 x 10“ 4 

Smooth mud flats 


O 

1 

Ln 

Sand flats 


4 x 10 -5 

Bare soils 


-4 

5 x io 

Thick grass tO.lO 

meters 

0.02 

Thin grass i 0 . 50 

meters 

0.05 

Thick grass <0.50 

meters 

0.09 

Field crops, tall 

grass ~ 1 meter 

0.10 - 0.14 

Forests, houses 


1.00 - 2.00 

Cities 


1.00 - 3.00 


i 


i 


Source: G. H. Fichtl, "Can Turbulent Drag Force 

Expansions Apart?" Paper presented at ASCE National Water 
Resource Conference Engineering Meeting, Atlanta, Georgia, 
January 24-28, 1972. 
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Continuity equation 


3u , 3v 3w 

3x 3y 3z 


0 


( 4 ) 


Momentum equations 


u 


3_u 

3x 


+ 




l dp 

P G 


+ 



fv 


(5) 


u 





+ 1 - 7 

P D 9y 9z yz 


fu 


(6) 


u 


9v 

3x 


+ 




^ 3 £ + 3 _ T 
P Q 3z 3z zz 



(7) 


where f is the Coriolis parameter, x is a shear stress, p Q 

and T q denote the density and temperature consistent with an 

adiabatic atmosphere. The stresses due to eddies lying in 

the horizontal plane (i.e., x , x ) have been neglected in 

xy yx 

comparison with the stresses arising from eddies lying in 

the vertical plane (i.e., x and x ). The bar over the 

xz y z 

dependent variables indicates an ensemble average. In con- 
sidering flow in the lower regions of the atmosphere certain 
assumptions and simplifications can be made in these equa- 
tions and will now be discussed in detail. 


Coriolis Effects in Atmospheric Shear Flows 

The terms in Equations 5 and 6 given by fu and fv 
represent the apparent forces produced by the Coriolis 
effect on atmospheric motions. Below 30-50 m the shear 
viscous forces predominate and the Coriolis forces are 
generally neglected. This is discussed by Tverskoi (15) 
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who argues that below a few tens of meters the integral term | 

I 

in I 

3 ! 

t(z) -t(z = 0) = f/ (u - u )dz i 

o 9 

\ 

obtained by integrating the equations of the geostrophic 
wind, is less than 10 per cent of the shear. For atmo- j 

spheric flow which occurs over much larger distances, the 
effects of the rotation of the earth must be considered 
however. The coordinate system on the earth is constantly 
accelerated in an inertial frame (fixed reference frame j 

relative to distance stars), hence, motions which are 
actually in straight lines in the fixed coordinate system 
appear from the ground to travel in a curved path. The 
fictitious Coriolis force in the relative frame of reference 

I 

accounts for this apparent deflection. In magnitude the 
Coriolis force is a function of the earth's angular velocity, 

Q, and the geographical latitude, 4> , and is proportional to 
the velocity of the motion considered. It deflects the 
velocity in a horizontal plane to the right in the Northern 
Hemisphere and in the opposite direction in the Southern 
Hemisphere acting perpendicular to the velocity vector, as 
shown by the' following relation: 

F c - [2fi sin cf>] V x k = f(V x k) 

where k is the unit vector in the vertical direction and f 
is the Coriolis parameter. 

With the magnitude of the Coriolis force determined 
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primarily by the velocity of the wind, its relative impor- 
tance in contributing to the momentum of the flow can be 
examined by considering the other dominant forces which act 
upon the atmosphere. Outside of the viscous region of the 
atmosphere, which generally extends to a height of about one 
kilometer, the motion of the wind is to a sufficient degree 
of approximation determined by a balance of the Coriolis 
force and the pressure gradient. However, as the effects of 
friction become significant within the boundary layer, the 
mean wind speed decreases causing the effects of the Coriolis 
force to also decrease. The balance of forces which deter- 
mines the momentum of the flow within the boundary layer 
includes the pressure field opposed by a combination of the 
viscous and Coriolis forces. Very near the surface, where 
the velocity approaches zero, the Coriolis force becomes 
negligible, and the atmospheric motion is determined by a 
balance of the pressure gradient and the viscous forces. As 
*a result of this change in the relative magnitude of the 
Coriolis and shear forces, two distinct regions can be 
identified within the atmospheric boundary layer: an inner 

layer extending from the ground to a height of about 50 
meters, in which the surface exerts the predominate influence 
on the flow; and an outer region whose structure is primarily 
determined by the Coriolis , pressure and shear forces . 

Therefore, in view of the physical situation which 
exists in the lower region of the atmospheric boundary layer, 
the Coriolis term in the governing equations for flow over a 
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surface obstruction can be neglected, permitting a simpli- 
fication of the equations to the following two-dimensional 
form: 


3u , 3w 
Sx 3z 


0 


( 8 ) 


u 


3u 

3x 


+ w 


3u 
3 z 


1 3t 

1 dp + xz 

p Q dx 3z 


vi 


3w 

3x 


+ w 


3w 1 3p 

3z 3z 

o 



(9) 

( 10 ) 


Approximation of the Reynolds' Stress by Turbulent Viscosity 
The turbulent boundary layer equations given by Equa- 
tions 8, 9 and 10 are obtained from the Navier-Stokes equa- 
tions by a method first proposed by Reynolds, in which he 
considered the instantaneous turbulent motion as a sum of a 
mean flow component and a randomly fluctuating component, 
given by 


u(t) = u + u' (t) , etc. 

Introducing the Reynolds' assumption into the Navier-Stokes 

equations and statistically averaging the instantaneous 

* 

motion over the ensemble of realization, the equations for 
general turbulent motion are obtained in which certain com- 
binations of the fluctuating component are contained as 
shown in the following relationships: 


3u 3w 
3x 3 z 


0 


* 


In practice the ergodic hypothesis is used and 
ensemble averages are estimated with time averages. 
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„ i- 3u _ 3p , . 3 2 u , 3 2 w 


-P. 


f o ' 


3x 


3u" r w " r 
32 


p o(“ i. + “Hi - 


■ H + 2u 


3 2 w 


3Z 


“ P 


3u'w ' 3w* 2 


►31 

i 

O 

3x 3z 

+ g 

j 


These equations are known as the Reynolds ' equations for 
turbulent flow in which the last group of terms on the right 
side can be interpreted as additional components of the 
stress tensor, called the Reynolds' stresses. Simplifying 
these equations, based on an order of magnitude analysis, 
results in the boundary layer equations given by Equations 
8, 9 and 10 in which the Reynolds' stresses are reduced to 
one term, pu'w', which is associated wxth an increase in the 
shear of the fluid. 

However,, in order to close this set of boundary layer 
equations, it is necessary to relate this additional stress 
to the mean flow. This is generally done by considering an 
analogy between the viscous shear which exists on a molecular 
scale and the turbulent shear resulting from the momentum 
transfer of macroscopic fluid particles. Accordingly, the 
Reynolds' stress is said to be the result of an eddy vis- 
cosity, e, which is related to the mean flow by an expres- 
sion analogous to the Stokes hypothesis for laminar shear 
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stress , 


T 


l 


pv 


du 

dz 


T = pu'w' 


pe 


du 

dz 


This approach in relating the turbulent fluctuations 
to the mean flow through the concept of the eddy viscosity 
has been used extensively in describing many turbulent flows 
however, it is limited in usefulness by dependence on 
empirical correlations of specific flow situations. 

Empirical Models of Eddy Viscosity 

The concept of the mixing length. The development of 
expressions for the eddy viscosity relating it to parameters 
of the flow, is generally based upon Prandtl's concept of 
the transport processes involved in turbulence by relating 
the turbulent fluctuations in velocity to the mean flow 
through the idea of a mixing length. 

Prandtl visualized turbulent flow as a random motion 
of macroscopic "lumps" of fluid which continually form and 
move through the fluid certain distances before being dis- 
persed into the surrounding fluid. The distance over which 
these fluid "lumps" retain their original properties before 
coming into equilibrium with their surroundings is termed 
the mixing length, l. A transfer of momentum, for example, 
can occur during this exchange process, resulting in a 
variation of the mean velocity at a particular location. 
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Consider two adjacent layers of fluid moving at different 
velocities, as shown in Figure 3, in which the mean flow is 
iji the x-direction. The mean velocity is specified as u at 
a distance z from a surface and can be approximated at the 
adjacent layer by the positive velocity gradient as u + it • 
(du)/(dz). Consider the motion of a fluid particle from the 
position z to z + l under the influence of a positive fluc- 
tuation- of velocity, w' , in the transverse direction. Since 
the momentum exchanged during this transfer of fluid results 
in a decrease in the velocity at z + l, Prandtl assumed that 
the momentum of the fluid particle remained constant so that 
the turbulent fluctuation in the hori zonal velocity was 
given by the difference in velocity of the two layers as 


u’ 


^ ( u + 




t 


du 

dz 


Further assuming that the transverse fluctuation, w' , is the 
same order of magnitude as u', it is given by 


w' = C • u ' 


Cfi. 


du 

dz 


Absorbing the constant C into the unknown mixing length, the 
Reynolds' stress is given by 


t = pu'w' 



du 

cfi" 


2 



du 

cU 


du 

dz 


( 11 ) 


with the eddy viscosity as 
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.2 du 
e = Si -3— 
dz 


In the neighborhood of a wall, the mixing length is 
generally assumed to be linearly related to the normal dis- 
tance from the surface 

l = KZ 

yielding the following expression for the eddy viscosity 


du 

2 2 du 

dz 

■ * z a? 


where k is the universal von Karman constant generally 
assumed as 0.4. Applying this approximation of the mixing 
length to the turbulent shear stress near the wall, gives 


2 2 du du 
T t = pK 2 d7 dF 


Assuming this shear stress is constant in a region next to 
the wall. Equation 13 can be integrated for the mean velocity 
profile 


u = — £n z + C 


The above logarithmic expression, as previously 
mentioned, is valid for turbulent flow over smooth surfaces, 
in which the mixing length approaches zero at the wall. 
However, for very rough surfaces the mixing length does not 
decay to zero near the ground, but approaches a size on the 
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and is given by the 


order of the roughness length, z q/ 
following relation 

£ = K (z + Z Q ) 


Therefore, the eddy viscosity over rough surfaces becomes 


e 


2 , ,2 

K (z + Z Q ) 


du 

37 


(15) 


leading to a logarithmic velocity distribution of the 
following form, previously introduced in Equation 3: 


u 



z + z 

o 


z 
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Since rough surfaces are characteristic of most natural 
terrain the eddy viscosity model used in the calculations 
carried out in this study is given by Equation 15. 

However, in verifying the numerical solution for 
turbulent boundary layer calculations, it is necessary to 
consider aerodynamic data of which a majority is obtained 
for flow over smooth surfaces requiring more sophisticated 
models of the eddy viscosity. Therefore, returning to Equa- 
tion 14, it is evident that very near the wall the predicted 
logarithmic velocity is not valid because it does not 
approach zero at the surface. This is a result of the fact 
that the dampening effect of the wall is not taken into 
account. In actuality there exists a very thin layer next 
to the surface in which turbulence is no longer a factor 


2 


I 


f 


and the shear stress is predominately determined by molec 
ular viscosity. In this region, known as the laminar sub 
layer, the eddy viscosity can be assumed zero. The thick 
ness of the sublayer is generally obtained from empirical 
relations (14) as approximately 


z 


sil 


= 5 



(16) 


Van Driest eddy viscosity model . The condition that 
the eddy viscosity abruptly becomes zero at a height above 
the surface given by Equation 16 is an oversimplification of 
the physical situation. To overcome the arbitrary nature of 
this assumption. Van Driest (17) proposed that the turbulent 
viscosity approach zero exponentially near the surface, by 
modifying the Prandtl mixing length as follows: 

i = kz[1 - exp(-z/A)] 


where 


/ V 7 p 

The Van Driest mixing length, valid for flow over a 
flat plate, was extended by Cebeci and Smith (18) to include 
the effects of a pressure gradient on the eddy viscosity. 
From the momentum equation, they approximate the shear 
stress close to the wall by 



21 


j 



Defining A as (26v) /(/t/p) , the mixing length can be 
expressed as 


t 

SL = kz j 1 - 
i 



d£ £ | 

dx p 




and the modified Van Driest model for eddy viscosity becomes 


e = 


2 2 / 

K Z \ 1 


- exp | - 


( T ^ 1/ 2 ] 

Z ' 1 

26v^p dx pj 


du 

dz 


(17) 


Although the above expression for the eddy viscosity 
is physically more realistic near the wall, in the outer 
regions of the boundary layer the surface no longer exerts 
the predominant influence on the flow. It is well docu- 
mented by experimental evidence (14) and (18) that at 
certain distances from the boundary, the velocity profile 
deviates from the logarithmic form given by the mixing 
length theory. This is a result primarily of a decrease in 
the influence of the wall on the characteristics of the 
turbulence, causing the outer region of the boundary layer 
to approach a free turbulent flow, in which the eddy vis- 
cosity becomes constant. Therefore, Clauser (19) assumed 
that in the outer region of the viscous flow the eddy vis- 
cosity is related to the overall parameters of the boundary 
layer by 


e 

o 


K U 
o e 


where 6* 


is the displacement thickness. 


(18) 
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However, due to the highly fluctuating' nature of the 
region separating the fully turbulent boundary layer from the 
relatively turbulent-free external flow, the constant vis- 
cosity given by Equation 18 should be modified by a damping 
factor experimentally determined from the intermittency 
factor of Klebanoff (20) as 


Y = 


1 + 5 - 5 (£j 


6 ’ 

) 


-1 


Separation of the inner and outer regions of the boundary 
layer is established by continuity of the eddy viscosity so 
that 


e 



if e • 
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e o 
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There are numerous other empirical models for eddy 
viscosity, but those considered in developing the turbulent 
boundary layer solution presented in this study are as 
follows : 

1. Prandtl mixing length model. 

2. Van Driest mixing length model. 

3. Van Driest model with the outer layer of Clauser. 

Conservation of eddy viscosity . A more recent theory 
for determining the eddy viscosity from mean flow parameters 
was advanced by Nee and Kovasznay (21) in which they con- 
sider the viscosity as a quantity which is conserved during 
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Lhe turbulent exchange processes. This method has been out- 
lined in Appendix C together with a numerical technique for 
incorporating it into the present boundary layer program. 
Preliminary investigation of this method has shown some 
promise; however, further work is necessary to determine its 
usefulness over the simpler mixing length models. 

Consideration of the Pressure Force in the Boundary Layer 
Equations 

Incorporating the assumptions introduced in the pre- 
vious sections concerning the Coriolis effect and the eddy 
viscosity, and applying the classical boundary layer order 
of magnitude analysis the equations become 


3u 3w 
dX Tz 


0 


( 19 ) 
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0 


( 20 ) 

( 21 ) 


These equations have been very successful in 
describing flow fields over aerodynamic surfaces, for which 
the viscous regions are very thin relative to the character- 
istic dimensions of the body. In fact, the derivation of 
the boundary layer equations from the Navier-Stokes rela- 
tions from an order of magnitude analysis requires the 
viscous region to be thin. However, despite the condition 
of a thin viscous region, boundary layer analyses have been 
applied to cases in which the boundary layer thickness was 
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on the order of the dimensions of the body but yielded good 
agreement with experiment, as shown in Reference (22) . This 
tends to indicate that the application of boundary layer 
theory is not limited to confined regions of a flow field 
near a body. 

One of the most significant results of the thin 
boundary layer assumption, as shown by Equation 21, is that 
the pressure field within the viscous region becomes inde- 
pendent of the coordinate normal to the flow. This result 
is important in establishing a solution of the thin boundary 
layer equations, because the pressure imposed on the flow 
can be obtained from the potential pressure distribution 
along the surface of the body and is a known quantity in the 
boundary layer equations . 

However, in considering atmospheric flow over a sur- 
face obstruction, the thickness of the boundary layer is 
generally many times greater than the height of the obstacle. 
Thus, the pressure disturbance produced by the body is con- 
tained within the viscous flow and no well defined outer 
inviscid flow field exists as in aerodynamic problems. This 
difficulty can be resolved, however, by considering the con- 
cept of the disturbed boundary layer introduced in Chapter I. 
Assuming that the region of influence which surrounds the 
body is contained within the atmospheric boundary layer, the 
equations which describe the motion over the obstruction are 
no longer dependent on the flow field which exists outside 
the viscous region, but become a function of the boundary 
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conditions. and pressure field which are given within the 
shear flow. 

Therefore, to adapt the boundary layer equations to 
very thick viscous regions containing variations in pressure 
normal to the flow, the streamwise pressure gradient of 
Equation 20 can be interpreted as some average pressure 
force which is representative of the distribution of this 
gradient through the boundary layer. In order to determine 
how the horizontal pressure field over the obstruction varies 
in the vertical direction and to examine the effects of this 
term on the equations, the potential pressure distribution 
along the streamlines over the elliptical cylinder will now 
be considered. 

Applying the potential flow theory outlined in 
Appendix A, the pressure and inviscid streamlines can be 
calculated as given in Figure 4 for the 2/1 aspect ratio 
ellipse. First, consider the pressure variation along the 
surface or zero streamline. Approaching from upstream, the 
pressure increases from the free stream value to a maximum 
along the leading edge and then rapidly decreases below 
ambient conditions to a minimum at the top of the ellipse. 
Since the potential flow is symmetric about the minor axis 
of the ellipse, the identical pressure distribution exists 
on the downstream side of the cylinder. However, the 
gradient of the pressure is not symmetric about this axis 
but is an odd function of the streamwise coordinate, x. 

This will become important in later discussions in 
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Figure 4. Potential pressure distribution over .the elliptical cylinder. 


considering the magnitude of the adverse pressure fields up- 
stream of the ellipse approaching the forward stagnation 
point and on the downstream face of the cylinder from the 
top of the ellipse to the trailing edge. It is obvious from 
Figure 4 that the latter pressure gradient is much more 
severe and will have a significant effect on the viscous 
flow. 

Applying the potential pressure along the surface to 
the boundary layer developed far upstream on the horizontal 
terrain, a qualitative picture of the viscous flow field 
over the body can be obtained. Approaching the obstruction 
from upstream, as illustrated in Figure 5, the large increase 
in pressure causes deceleration of the flow and thickening 
of the internal boundary layer. As the pressure gradient 
becomes increasingly more severe, the momentum of the fluid 
becomes very small until the boundary layer is no longer able 
to remain attached to the surface, (point A) . Further down- 
stream, the flow is accelerated by the decrease in pressure 
over the front face of the ellipse causing the boundary 
layer to reattach to the body at point B. From this point 
the flow continues to accelerate to the top of the cylinder 
where the adverse pressure field on the rear side again 
forces the boundary layer to separate at point C. Much 
further downstream the pressure field begins to recover from 
the severe gradients discussed above and returns to ambient 
conditions causing the flow to reattach once more (point D) 
and readjust to the conditions along the horizontal surface. 
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Figure 5. Qualitative description of the flow field over the elliptical cylinder 



The areas of the flow field enclosed by streamlines AB up- 
stream of the ellipse and CD on the downstream side are 
separation regions in which the flow is recirculated, 
tending to become stagnant with the velocity near the sur- 
face i,n the opposite direction to that in the outer flow as 
illustrated by the profiles in Figure 5. This separated 
flow region is an extremely complicated motion which cannot 
be predicted by a boundary layer analysis. Therefore, 
additional assumptions must be made about the adverse pres- 
sure field which produces this separated flow to enable the 
boundary layer solutions to be computed over this region to 
the top of the ellipse. 

Comparing the variation of the streamwise pressure 
distribution along different streamlines above the surface, 
as shown in Figure 4, page 27, it can be observed that as 
the distance from the body increases, the severity of the 
pressure gradient decays until outside a certain streamline, 
the adverse pressure upstream of the body no longer causes 
separation of the flow. Since the pressure field near the 
surface primarily determines the flow within the separation 
region, it is reasonable to assume that the pressure gradient 
given along this nonseparating streamline is more character- 
istic of the overall flow which actually negotiates the 
obstruction. With this in mind, the boundary layer analysis 
proposed in this study was carried out with the assumption 
that the flow outside the separation region is governed by 
an average pressure field given by at calculated along the 
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first nonseparating streamline over the elliptical body. 

Improved Approximation of the Vertical Pressure Field 

The assumption of an average pressure gradient which 
is constant in the direction normal to the flow is expected 
to be a good approximation of the conditions which exist in 
a region near the body, as will be substantiated in later 
discussions concerning the results of calculations carried 
out over the elliptical cylinder using this mean pressure 
distribution. However, if interest lies in flow fields much 
further removed from the surface obstruction, a more physi- 
cally realistic pressure field is suggested by the fact that 
the pressure gradient approaches zero far from the body. 
Therefore, for boundary layers which extend well above the 
obstruction, an improved pressure field would be one which 
decays in the vertical direction from the initial value 
determined along the nonseparating streamline to zero as the 
distance from the body approaches infinity. This type of 
pressure variation will permit a matching of the flow within 
the region of influence of the body to the outer undisturbed 
logarithmic profile. Under this assumption, the form of the 
pressure distribution introduced into the boundary layer 
equations is given by 


dp (x, z) 
dx 



• q[ z) 


where [(dp)/(dx)] is the pressure gradient given along the 

ns 

nonseparating streamline, and q(z) is the vertical decay 
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function, which will be discussed in more detail in a later 
section . 

Boundary Conditions 

The governing equations described in the preceding 
sections are complete only when the proper boundary condi- 
tions along the surface and matching conditions with the 

» 

outer flow are specified. 

Lower boundary conditions . Along the lower boundary 
the no-slip condition will generally apply, and thus the 
mean velocities and turbulent fluctuations approach zero 
near the ground. However, the roughness associated with the 
surface also has an influence upon the flow in the vicinity 
of the lower boundary. One effect is to increase the tur- 
bulent motion, as noted previously, by assuming a larger 
mixing Length at the surface. An additional effect of the 
rough surface is caused by the separation of the flow on a 
very small scale from the individual roughness elements 
which increases the dissipation of momentum due to larger 
form drag over the surface. Since neither of these effects 
are accounted for by the simple no-slip condition, it is 
assumed that the logarithmic velocity distribution given by 
Equation 3 applies in a very thin layer next to the surface 
in which the shear stress is assumed constant, as shown in 
Figure 6. With this boundary condition the influence of the 
surface roughness on the velocity can be incorporated into 
the solution of the boundary layer equations. 
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Outer boundary conditions . The outer boundary condi- 
tion usually imposed on classical boundary layers is the 
potential velocity along the surface streamline of the body. 
This assumption is based on the fact that the boundary layer 
is extremely thin and the surface curvature is not large. 

For the problem at hand, however, the pressure distri- 
bution over the body is not assumed to be that corresponding 
to the potential flow along the surface streamline but that 
corresponding to the first streamline for which separation 
at the front of the ellipse does not occur. The outer 
boundary condition is therefore taken as the potential 
velocity, U e , along this first non-separating streamline. 
Physically one defends this assumption by arguing that the 

flow which negotiates the body avoiding the separation region 

dU. 

is driven by some average pressure characterized by 


U 


e dx J ns 

and the velocity scale for the disturbed shear region or 
internal boundary layer is Ue ng • Although it is known that 
far from the body the velocity must return to the logarithmic 
profile of the undisturbed flow it is expected that Ue ns may 
represent the peak in the velocity profile and the solution 
corresponding to this boundary condition might be valid very 
near to the wall. Later refined calculation shows this 
region extends to a height on the order of z/b = 0.5. A 
number of solutions were carried out utilizing this boundary 
condition and pressure distribution and are expected to be 
representative of the flow very near the body. 

However, once again considering flow fields at 
larger heights, the outer boundary condition must be 
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reconsidered in light of the rapid decay of the pressure 
field which exists in the greater region of influence above 
the surface obstruction. Since the flow tends to return to 
the logarithmic form as the effects of the obstacle become 
negligible, a more realistic outer boundary condition 
matching the outer undisturbed flow -with the boundary layer 
developed over the body would be given by the logarithmic 
velocity distribution. Therefore, in considering the flow 
field in a larger region above the elliptical cylinder, the 
outer boundary condition used in connection with the decay- 
ing pressure field is given by the logarithmic velocity 
characteristic of the undisturbed flow at an arbitrarily 
selected elevation. 

The Effect of Separation Regions on the Boundary Layer 

The separation regions which exist upstream and down- 
stream of the elliptical cylinder, as shown in Figure 5, 
page 29, are due to the strong adverse pressure gradients 
near the surface which have been neglected in the initial 
solutions to be presented because the pressure field is 
approximated along the nonseparating streamline. For this 
reason the boundary layer equations can be integrated over 
the forward separation region without regarding the effects 
it produces on the flow. However, this method cannot be 
extended over the rear separation bubble because the adverse 
pressure gradient along the same streamline becomes much 
more severe on the downstream side of the body and prevents 
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the compul al innal scheme from negotiating Lins reg ion. As 
discussed previously in comparing the magnitude of the 
positive pressure gradients shown in Figure 4, page 27, the 
much larger pressure forces occurring downstream signifi- 
cantly affect the boundary layer flow, so that in order to 
apply the present integration technique over the rear 
separation region, the pressure field must be obtained from 
an in viscid streamline very far from the body. Since the 
pressure gradients downstream will be greatly reduced, the 
solutions are no longer believed to give a meaningful 
description of the disturbed flow field created by the sur- 
face obstruction. Therefore, the solutions given in this 
study for flow over the elliptical body are carried out only 
to the top of the cylinder. 

To incorporate some of the influence of the separa- 
tion regions on the pressure field near the surface and to 
obtain a more realistic boundary condition on the velocity 
over these regions, some preliminary solutions have been 
computed with the recirculating flow of the separation 
bubbles simulated by stationary, inviscid vortices calcu- 
lated from potential flow theory (23), Although not dealing 
directly with the flow within these regions, this technique, 
outlined in Appendix B, gives an improved approximation of 
the pressure field imposed on the flow by the elliptical 
cylinder. A typical result of this calculation is given in 
Figure 7, showing the inviscid streamlines over a 2/1 aspect 
ratio ellipse with attached vortices together with the 
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potential velocity distribution along these streamlines. in 
addition, the potential velocity along the nonseparating 
streamline over the same body without vortices is also shown 
for comparison purposes. 

Despite the improved approximation of the pressure 
field obtained with this model which includes the effects of 
the separation regions, the adverse pressure gradients 
associated with the flow approaching the stagnation points A 
and B are still sufficiently severe near the body to prevent 
the boundary layer calculation from negotiating these 
regions. Therefore, the applicable pressure distribution is 
again assumed to be given along an outer streamline over the 
vortices for which the pressure gradient does not produce a 
condition of separation in the boundary layer solution. 
Although this pressure field still does not simulate that 
given along the surface, it does include a part of the 
significant influence of the separation regions on the 
boundary layer which is neglected by the initial approxima- 
tion of the pressure gradient over the elliptical cylinder. 

In addition to the effect of the vortices on the 
pressure field, the potential velocity along the streamline 
dividing the separation region from the outer boundary layer 
can be incorporated into the flow field developed over this 
region as a lower boundary condition. Since the velocity 
which physically exists in the vicinity of the dividing 
streamline does not approach zero, the slip-velocity provided 
by this potential solution is a more realistic measure of 
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the lower boundary condition affecting the equations than 
that previously introduced for the calculations made over 
this region. However, the potential velocity calculated 
over the vortices, as shown in Figure 7 , page 37/ approaches 
free stream values at its maximum between the stagnation 
points on the dividing streamline. Since this characteristic 
of the potential solution is a result of the effect of the 
vortex on an inviscid fluid, the viscous dissipation which 
actually occurs in this region decreases the potential 
velocity, and therefore, the lower boundary condition 
assumed from this model will be a fraction of the inviscid 
velocity, which has been chosen for the calculations made 
with this method as 0.25. 

Curvilinear Effects of the Streamline-Oriented Coordinate 
System 

Since the pressure gradient imposed on the boundary 
layer is given along the inviscid streamlines over the 
elliptical cylinder, the coordinate system used must also be 
oriented along these streamlines resulting in a curvilinear, 
orthogonal system, as shown in Figure 8a. The curvature of 
this coordinate system is small throughout the region sur- 
rounding the ellipse, except in the vicinity of the upstream 
and downstream stagnation points where the slope of the sur- 
face is discontinuous. This is important in view of the 
fact that the boundary layer concept is generally assumed to 
be a first order approximation of viscous flow and neglects 
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Figure 8. Coordina 


the higher order effects produced by body curvature. There- 
fore, the accuracy of the present boundary layer analysis is 
expected to deteriorate in the vicinity of the stagnation 
regions but be relatively good upstream of these regions and 
on top of the ellipse. 

Moreover, since the flow field within these stagna- 
tion regions has been neglected by the assumption, pre- 
viously discussed in Chapter II, of a pressure distribution - 
imposed on the flow as given by that along the nonseparating 
streamline, there is little value in modifying the boundary 
layer equations with the additional complications of a 
curvilinear coordinate system. Consequently, the calcula- 
tions have been carried out in Cartesian coordinates (x,y) 
as shown in Figure 8b. The x-direction is measured along 
the inviscid streamline with the y-axis extending perpen- 
dicularly to it at each x-point. The velocity profiles 
calculated in this coordinate system are then assumed to 
exist physically on the surface, which is a reasonable 
assumption because the geometrical effects of the body on 
the flow enter the equations only through the pressure 
gradient. However, because of the assumptions made about 
the pressure distribution and curvature effects over the 
stagnation regions, the interpretation of the velocity pro- 
files calculated in these areas of the flow field is diffi- 
cult to apply physically. Consequently, the profiles 
presented within this investigation are those in the regions 
of flow upstream and on top of the ellipse where the 
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solutions are considered to be the most meaningful, 
neglecting the profiles calculated over the separation 
regions . 
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CHAPTER III 


NUMERICAL SOLUTION OF THE TURBULENT 
BOUNDARY LAYER EQUATIONS 


The turbulent boundary layer equations together with 
the boundary and initial conditions and the assumption of an 
eddy viscosity model to relate the turbulent motion to the 
mean flow variables, form a closed set of nonlinear, para- 
bolic, partial differential equations given by Equations 19 
and 20. This set of equations must be solved by numerical 
methods which will now be introduced. Approximating the 
pressure gradient term by the inviscid pressure distribution 
along the nonseparating streamline, the Bernoulli equation 
can be used to relate this term to the potential velocity as 
follows : 


1 dp 
p dx 


U 


dUe 
e dx 


( 22 ) 


Adopting some characteristic length, L, and velocity, U^, 
from the flow field over the surface obstruction, these 
equations can be nondimensionalized into the following form: 
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where the dimensionless quantities arfe indicated by stars as 
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and Re denotes the dimensionless parameter defined as the 
Reynolds 1 number by 
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oo 


V 


Dropping the dimensionless notation and expanding, the 
equations to be solved are given by 
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A transformation of the vertical coordinate, z, which 
will compress the vertical length scale near the wall where 
the velocity gradients are the most severe, is now intro- 
duced into the governing equations in the form of a stretch- 
ing function c - f(z). Applying this coordinate transforma- 
tion to the vertical derivatives results in the following 
relations : 


3 _ 

3z 


f ' 



44 



— — y « ff ’ (z) ] 2 + f"(z)^- 

9z Z 9C 3C 


Defining the functions 

* x (0 = f ’ (z) = f'tg(c)]; 4 2 (c) * f h [g(C)] 


the boundary layer equations become 
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The stretching function, £(z), used in the calcula- 
tions presented in this investigation, is given by the 
following quadratic relation: 

2 

z = Ic, + me + n 


in which the constants £ , m and n were adjusted to most 
efficiently space the calculation points through the boundary 
layer. 

The Implicit, Finite Difference Technique for Solution of 
the Boundary Layer Equations 

Assuming that the flow field over the elliptical 
cylinder is represented by the Cartesian grid given in 
Figure 86, page 40, the equations of motion cam be solved by 
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approximating derivatives by finite differences expressed 
in terms of the nodal points of this grid. This numerical 
approximation results in a set of finite difference equa- 
tions involving the unknown velocities at nodal points along 
an x-station, m + 1, in terms of known values at the pre- 
ceding station m. Specifically, the equations are written 
in finite difference form at each point of the unknown m + 1 
station using the following expressions for the x- and z- 
deri va fives : 
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However, since the momentum equation of the boundary 
layer is nonlinear in the inertial terms, an iterative pro- 
cedure must be introduced into the above numerical method. 
Initially, the values of the velocity at the known station m 
will be used to linearize the equations at m + 1. Numeri- 
cally solving the system of equations at m + 1, this solu- 
tion can then be used to linearize the momentum equation for 
a. second iteration. Continuing this process by repeatedly 
linearizing the boundary layer equation with successive 
solutions for the velocity at m+1 until convergence is 
achieved, results in the numerical solution of the flow 
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field at this x-station. Proceeding downstream to the next 
station, m + 2, this method can be repeated by using the 
converged solution at m + 1 to initially linearize the equa- 
tions at m + 2 continuing the above iteration procedure 
until convergence of the solution at m + 2 is reached. In a 
similar manner this numerical procedure can be marched down- 
stream over the entire flow field obtaining the turbulent 
boundary layer solution. 

Returning to the details of the numerical solution, 
the iterative technique and the finite difference approxima- 
tions can be introduced into Equation 28, resulting in the 
following algebraic expression of the boundary layer equa- 
tions at a point (m + 1, n) : 
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where the superscript T indicates a trial value from the 
preceding iteration. Rearranging into the following general 
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form. Equation 29 becomes 


fl n u n-l + Vn + Wl = V 2 - " ‘ 1 


where the expressions A , B , C and D are given by 

n n n n J 


w T , • * . *2L_ 

m+ 1 , n 1 2 A C 


1 + e 


m+1 , n 


+ cl 

.2 ^2 


Ax_j 

2AC ! 


Re 13C j 


m+1 ,n 


u , . + — 

m+1 ,n 


1 + e 1 . 

m+l,n . 2 

Re V 1 


-w , , • 4> 

m+ 1 , n 1 


1 + e 


m+l,n 


£x_| + _1 

2 A t; Re 


m+ 1 , n 


D - u , . * u + U ———[ * Ax 

n m+1 ,n mn e dx , , , 

'm+1 


Applying these equations at each of the nodal points 
between the surface and outer boundary conditions, results 
in a system of N - 2 algebraic equations as mentioned above. 
This set of simultaneous equations yields a tridiagonal 
coefficient matrix which can be solved by an efficient 
matrix inversion technique involving the following recursion 
formula: 


u n = .i + F or 

n n n+1 n 


u . =? E . u + F . 
n-1 n-1 n n-1 
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Substituting into Equation 30, 


A u , , + B u +C[E .u + F ,]=D 
n n+1 n n n n-1 n n-1 n 


and solving for u 


n 


-A u , , D - C F , 

n n+1 n . n n-1 

n B+CE, B+CE, 
n n n-1 n n n-1 


the unknown coefficients E and F become 

n n 


-A 


E = 


n 


n B + C E . 
n n n- 1 


F = 


D - C F . 
n n n-1 


n B + C E . 
n n n-1 


(32) 


(33) 


which with Equation 30 define the velocity profile at sta- 
tion m + 1, providing the proper boundary conditions are 
specified. 

The lower boundary condition of a no-slip velocity at 

the wall is incorporated into Equation 31 to yield values of 

E. and F, of zero. Using these initial values E and F can 
11 3 n n 

be calculated through the boundary layer by Equations 32 and 
33. The outer boundary condition is used in Equation 31 to 
initiate the calculation of the velocities toward the wall- 
However, in order to incorporate the conditions imposed on 
the velocity by the surface roughness, the above numerical 
procedure must be altered in calculating the velocities near 
the wall by applying the logarithmic velocity distribution 
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to the points closest to the lower boundary. Recalling that 
the logarithmic law implies a layer of constant shear, as 
shown in Figure 6, page 33, the friction velocity, u* , will 
also be constant in this region. Therefore, assuming the 
value of u* is given by 


u* 


<u ( 3) 

fz(3) ♦ Z Q| 


the velocity at point 2 can be calculated by the logarithmic 
law as 


u. 


u { 2 ) 


In 


[ z ( 2 ) + z i 
o I 


) 


and u(l) will remain zero at the wall. 

In summary, the method for calculating the velocity 
profile at m + I is: 

1 . Calculate A , B , C and D from known values, 
n' n n n 

2. Calculate E and F by Equations 31 and 32 

n n ^ 

starting with n = 1 to n = N - 1, with = 0, 

*i ^ °- 

3. Calculate u(n) inward from the outer boundary 
condition to n = 3 by Equation 30. 

4. Calculate u(2) by the logarithmic law. 

The velocity profile in the vertical direction can be 
obtained by integrating the continuity equation given by 
Equation 27 using the velocity profile calculated above. 
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Approximating this integral by the trapezoidal rule, the 
vertical velocity component can be calculated by successively 
solving the following equation from the wall to the outer 
edge. 


w , , ( c ) 
m+1 ^n 




ll 1 


i' u m+l,n-l u m,n-lj 
Ax 1 


+ 


f u m+l,n U m,nj 

! Z* I 


Az 


The above equation together with Equations 30, 31 and 
32 and the boundary conditions give the numerical solution 
of the boundary layer at station m+1, after a suitable 
model of the eddy viscosity and pressure distribution have 
been introduced. 


Comments on Convergence and Accuracy o f the Nume r ical 
Solution 

In developing the numerical solution of the turbulent 
boundary layer equations, outlined in the previous section a 
balance between the convergence and accuracy required in the 
solution and the amount of computing time necessary to meet 
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these conditions must be consider -d. One of the most 
important factors affecting this balance is the grid spac- 
ing, Ax and Az, employed in approximating the derivatives by 
finite differences. Therefore, the following section will 
discuss the effects of this numerical increment size on the 
convergence, accuracy and computing time and outline some of 
the developmental work done in optimizing the step size. 

Convergence of the iterat i ve procedure . The iterative 
procedure, previously introduced to solve the numerical 
equations, must converge in order to obtain a meaningful 
solution of the boundary layer; that is, the difference 
between successive approximations of the solution at a 
particular x-station, m + 1, must continually decrease as 
the number of iterations increases. However, the convergence 
of the solution to the exact value, in which the difference 
between successive calculations is zero, requires a very 
large number ot iterations increasing the costs in computing 
time well beyond the value of the improvements made in the 
solution. Therefore, it is necessary to introduce a small 
degree of inaccuracy due to incomplete convergence as a 
compromise for. reasonable computing time. Thus the optimi- 
zation of a numerical scheme involves reducing the number of 
iterations f increasing the rate of convergence! required to 
efficiently obtain a given level of accuracy. 

In any iterative technique, the accuracy of the 
initial guess made at the converged solution determines the 
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number of iterations required to reach that solution. In 
the present method the initial guess, used in linearizing 
the equations, is the calculated velocity profile at the 
preceding x-station. Therefore, the strength of the veloc- 
ity gradient in the x-direction significantly influences the 
rate of convergence and determines the size of the finite 
increment in the x-coordinate, as shown in the next section. 

The effect of truncation error on accuracy . In addi- 
tion, to the errors brought about by incomplete convergence, 
the truncation error, or the difference between the solution 
of the finite difference equations and that of the differ- 
ential equations, is also a problem associated with a 
numerical solution. Since the derivatives are approximated 
by incomplete series expansions of differences made over 
finite distances, the truncation error is a function of the 
grid spacing, and ffects not only the rate of convergence 
but also the accuracy of the solution, becoming most critical 
in regions where the gradients of velocity are the largest. 
With this in mind, the initial investigation of the turbulent 
solution indicates that the maximum horizontal gradient 
occurs in regions where changes in pressure are the most 
severe, that is, near the separation regions upstream and 
downstream of the obstruction, and that the maximum velocity 
gradients in the vertical direction occur near the wall. 

Optimization of the numerical scheme in approximating 
horizontal gradients in regions of large pressure variation 
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has been achieved by allowing the horizontal increment 
spacing, Ax, to vary as a function of the number of itera- 
tions required for convergence at a particular x-station. 

As the iterations increase. Ax is decreased so that the 
number of iterations is kept below an upper limit. However, 
for large magnitudes of the adverse pressure gradients, 
which occur near the separation regions, the iterative pro- 
cedure will not converge, so that control on Ax will 
continue to decrease the step size, effectively determining 
a separation point of the boundary layer. 

The truncation errors associated with the approxima- 
tion of the vertical velocity gradients are also controlled 
by a variable grid spacing in the vertical coordinate, 
introduced previously by the stretching function, t(z). 

Since the velocity gradients are largest at the wall, the 
most important factor in determining the magnitude of the 
truncation errors js the physical grid size at the wall, 
denoted by (Az) ^ In the present numerical scheme this 

value is constant throughout the calculation field and is 
determined as a fraction of the boundary layer thickness of 
the velocity profile used as an initial condition on the 
equations . 

The following section presents a portion of the work 
carried out to optimize the numerical technique through an 
efficient choice of step sizes. Toward this end numerical 
solutions are compared with experimental data for turbulent 
flow over a flat plate. 
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Comparison of Numerical Solutions with Known Turbulent Flows 

Comparison of the turbulent boundary layer solution 
with the measured mean velocity profiles of Wieghardt (25) 
for flow over a smooth flat plate without imposed pressure 
gradients is given in Figure 9 , showing several calculated 
velocity distributions at one x-location for different 
values of the vertical increment size A£ . All calculations 
were started at an initial station, x^, well upstream, using 
a profile interpolated from the data at this point, and were 
extended downstream using a constant Ax step size of one 
centimeter, utilizing the Van Driest-Clauser eddy viscosity 
model, previously discussed. 

From earlier discussions, the critical factor in 
limiting the truncation errors associated with the vertical 
velocity gradient was determined to be the physical step 
size at the wall, (Az) z _q, which in all calculations is set 
by the stretching function as 10 per cent of A q . The 
importance of this factor in determining the accuracy of the 
solution is dramatically illustrated in Figure 9 by the 
obviously erroneous velocity profile obtained in the calcu- 
lation employing (Az) z _q of 0.2 millimeters, which represents 
approximately 2.0 per cent of the initial boundary layer 
thickness. The other two profiles in this figure, calcu- 
lated with smaller values of (Az) A of 0.5 per cent and 1.0 

z=0 

per cent of the initial boundary layer thickness, respec- 
tively, are well behaved with the smallest value giving a 
slightly better approximation of the data. Figure 10 
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on the calcula- 


illustrates the same influence of (Az) _ 

z=u 

tions of the displacement thickness over the plate. 

Based on the above discussion, it was decided to 
limit the maximum value of (Az) z _q to 1.0 per cent of the 
initial boundary layer thickness, implying an upper limit on 
the numerical step size, AC, of 10 times this value. 

Optimizing the rate of convergence, as previously 
discussed, by controlling the number of iterations through 
variation of Ax, basically determines the size of this 
increment. However, Ax was not allowed to grow to a magni- 
tude of larger than 10 times AC in order to control the 
truncation error introduced by the approximation of the 
horizontal gradients. 
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CHAPTER IV 


DISCUSSION OF RESULTS FOR THE FLOW FIELD 
NEAR ELLIPTICAL OBSTRUCTIONS 

Numerical solutions of the turbulent boundary layer 
equations have been carried out for* a number of different 
cases involving atmospheric flow over elliptical cylinders 
of various aspect ratios and surface roughnesses . The 
results of these calculations are divided into two types: 
the first, presented in this chapter, involves the analysis 
of the flow field in a region above the obstruction on the 
order of one ellipse height in which the pressure gradient 
and outer boundary condition are specified along the non- 
separating streamline; and the second, to be presented in 
Chapter V, in which the analysis is extended to include the 
flow field in a larger region of influence above the 
elliptical cylinder, incorporates the decaying pressure 
distribution and outer boundary conditions dictated by the 
logarithmic velocity - 

Considering the flow field close to the body, the 
following solutions are obtained from a region of calcula- 
tion shown in Figure 11. The calculation procedure, 
described in Chapter III, is initiated with the logarithmic 
velocity profile, characteristic of the approaching wind, at 
a position upstream of the body of 20 times the height of 
the ellipse. This location is chosen because it represents 
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Figure 11. Region of solution for flow over an elliptical cylinder. 


a distance outside of which the pressure disturbance created 
by the obstruction is less than one per cent of the free 
stream pressure. Proceeding downstream, the velocity pro- 
files are calculated as shown, using the potential velocity, 
U fi , as an outer boundary condition, passing through the 
region of adverse pressure where the velocities decrease and 
the boundary layer grows to the favorable pressure on top 
of the ellipse where the flow is accelerated. 

The results of this calculation of the boundary layer 
over the elliptical obstruction can be conveniently dis- 
cussed in terms of three relevant parameters which evolve 
from the dimensionless governing equations. Inspection of 
the boundary layer Equations 25 and 26 and the assumptions 
concerning the pressure field and boundary conditions which 
apply, identifies these parameters as: first, the aspect 

ratio, k, of the elliptical cylinder, defined as the ratio 
o.f the length major axis which is parallel to the ground, to 
the height of the ellipse which has been fixed in all cases 
considered to a value of one meter; second, the surface 
roughness characterized by a vertical length scale, z q ; and 
third, the Reynolds number, expressed in terms of the height 
of the elliptical cylinder, b, and a characteristic velocity, 
U w , defined at this same height by the logarithmic velocity 
law. 

The Effect of Elliptical Aspect Ratio on the Boundary Layer 

The effects of a parametric variation of the aspect 
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ratio on the velocity profiles at the top of the ellipse are 
shown in Figures 12 and 13 for two different surface rough- 
nesses. The apparent unusual behavior of the curves at the 
coordinate position (1, 1) can be explained both by the 
nature of the boundary condition chosen and by the smoothing 
effect of the numerical integration. In the former case 
u/U e becomes unity at an elevation slightly greater than 
z/b = 1 due to the displacement of the potential solution 
streamlines along which U e is evaluated. The forcing of the 
solution by the boundary conditions makes the solution in 
this region less reliable than the remaining portion of the 
velocity profile. In both cases, the effect of increasing 
the aspect ratio is a decrease in the mean wind velocity at 
the top of the ellipse. It is also of interest to compare 
the changes produced by the different ellipses in the 
initial logarithmic velocity profile existing far upstream. 

For the low aspect ratios, the velocity profiles show wind 
speeds at the top which are greater than the approaching 
wind, but the larger aspect ratio produces velocities which 
are very close to the initial logarithmic distribution, 
characteristic of the undisturbed flow, deviating signifi- 
cantly only in the outer portion of the profile where the 
boundary layer growth causes a decrease in velocity from the 
initial values. 

The effect of the elliptical aspect ratio on the veloc- 
ity profiles can be explained in terms of the pressure field 
assumed to exist over the cylinders as shown in Figure 14. 
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Figure 12. Effect of elliptical aspect ratio on velocity 
profile at the top of the ellipse for z Q /b = 0.005. 
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re 14. Comparison of the pressure distribution over 
elliptical cylinder. 


Comparing the adverse pressure gradient given along the 
nonseparating streamlines for the large and small aspect 

ratio, it is apparent that the pressure distributions are 
very similar with the 2/1 ellipse having only a slightly 
larger gradient. On the other hand, the favorable pressure 
gradient is obviously quite different for the two cases. 

For the 10/1 ellipse the gradient is large only for a short 
distance near the leading edge then rapidly decreases to a 
very small value over a significant portion near the top of 
the ellipse. For the 2/1 ellipse the favorable gradient is 
much larger approaching zero only in a very small region at 
the top of the cylinder. The additional acceleration given 
to the flow by this larger pressure force is the primary 
cause of the higher wind speeds shown in Figures 12 and 13 
at top of the smaller aspect ratio ellipses. In addition, 
the fact that the velocity tends to return to the logarithmic 
distribution at the top of the longer ellipse ca p. also be 
explained by considering the region of a favorable pressure 
gradient. Recalling that the logarithmic velocity is 
characteristic of zero pressure gradient flow, it is not 
surprising that the boundary layer at the top of the 10/1 
ellipse approaches the logarithmic form of the upstream, 
undisturbed wind profile. 

The £ffe c t of Sur f ac e Roughness on the Boundary Layer 

The effects of variation in surface roughness on the 
velocity of the wind at the top of the 2/1 and 10/1 ellipses 
are shown m Figures 15 and 16, respectively, for several 
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Figure 16. Effect of surface roughness on the velocity 
profiles at the top of 10/1 ellipse. 
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different roughness values, z /b, ranging from 0.005 to 
0 020. The decrease in the velocity at the top of the two 
ellipses caused by increasing the roughness, can readily be 
explained by examining the influence of the individual 
roughness elements on the flow in the vicinity of the ground. 
As discussed previously, there are two effects produced by 
the protrusions of the surface into the flow: the first is 

the increased drag of the fluid on the surface caused by the 
separation of the boundary layer on a small scale; and the 

second is the larger shear brought about by the more turbu- 
lent eddying motion produced in the region close to the 
surface. Both of these effects are enhanced by larger 
roughness creating a greater dissipation of the momentum 
through the boundary layer. 

Comparing the effects of an increase in roughness on 
the two different ellipses shown in Figures 15 and 16, it 
appears that the roughness change does not influence the 
flow fields for the two different aspect ratios to the same 
degree. This is more apparent in Figure 17 where the 
difference, Au, in velocity calculated with roughness z Q /b = 
0.005 and that for z Q /b = .020 is plotted as a function of 
height above the surface. From this figure it is obvious 
that increased roughness has more influence on the boundary 
layer near the surface of the 2/1 ellipse than the 10/1, but 
that this influence decays rapidly toward the outside of the 
surface layer causing roughness effects to become more 
influential in the outer region for the 10/1 ellipse. 
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These results can be explained by again comparing the 
pressure field over the top of the 2/1 and 10/1 elliptical 
cylinders. As previously discussed the pressure gradients 
which exist along the top of the larger aspect ratio 
ellipse are extremely small, so that the force created by 
this gradient is negligible, and the momentum of the flow in 
regions near the top of the 10/1 ellipse is effectively 
determined by the viscous forces only. Therefore, the 
viscous effects produced by the surface roughness extend 
much further into the boundary layer for the 10/1 ellipse, 
explaining the larger values of Au calculated in the outer 
region of the boundary layer as observed in Figure 17, 

On the other hand, the much larger favorable pressure 
gradients which exist over the top of the 2/1 ellipse have 
an influence on the flow which can be used to explain the 
effects of a change in roughness observed near the surface. 
Because of the larger acceleration produced by these favor- 
able pressure forces, the velocity gradient near the wall is 
greater for the 2/1 ellipse than for the 10/1 ellipse. 
Therefore, since the viscous shear is directly related to 
this gradient, the effects of roughness will also be larger 
in this region, leading to the results shown in Figure 17. 

The Reynolds 1 Number as a Parameter 

Nondimensionaliz ing the governing equations by a 
characteristic length and velocity produces the Reynolds' 
number as a parameter of the boundary layer. However, by 

defining the Reynolds' number in terms of the logarithmic 
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wind profile which exists far upstream, it. is no longer an 

independent parameter of the equations but becomes a func- 

/ 

tion of the surface roughness, as shown in the following 
equation: 


Re 


u » b 

v 



In 


i'b + z 


o 



(34) 


On the other hand, this definition of the Reynolds’ number 
introduces a new variable into the problem as the initial 
friction velocity of the undisturbed flow, u* Q , which has 
been chosen as the independent parameter, entering the 
equations through the Reynolds’ number. 

The influence of the initial friction velocity on the 
overall flow field is illustrated in Figure 18 by the dis- 
tribution of the displacement thickness over the 2/1 aspect 
ratio ellipse. Two curves are shown, one for a value of 
u* Q = 0.40 meters/second and the other for u* Q = 1.00 
meters/second which correspond to Reynolds' numbers of 
3.8 x 10~* and 9.5 * 10^, respectively. For the values of 
friction velocity used in these two calculations, the 
maximum change occurring in the displacement thickness for 
an increase in u* by a factor of 2.5 is only 0.87 per cent. 
Thus, it is suggested that the Reynolds' number, on the order 
of the magnitudes given above, has little effect on the 
boundary layer. This conclusion can also be verified by 
examining the nondimens lonal boundary layer equation 
(Equation 24). Expressing the eddy viscosity derived from 
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the mixing length concept as given in Equation 15 in non- 
dimensional terms. 
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and substituting into Equation 24 results in the following 
expression for the boundary layer equation: 
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It can be seen from this equation that the effect of the 
Reynolds' number will be negligible through a majority of 
the boundary layer if the Reynolds' number is large as is 
characteristic of atmospheric flows. 

Although the overall effect of the Reynolds' number 
is small, il us also evident from Equation 35 that very near 
the wall where the turbulent shear is decreased, the 
importance of the Reynolds' number is increased. This is 
clearly demonstrated in Figure 19 by comparing the changes 
in the velocity profiles produced by the Reynolds' numbers 
considered above, A velocity difference, u - u', is plotted 
as a function of height above the surface, where u is the 
wind velocity calculated with an initial friction velocity 
of 1.00 meters/second, and u' is the wind velocity for a 
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Figure 19. Effect of the initial friction velocity and Reynolds' number on the 
displacement thickness. 




friction velocity of 0.40 meters /second. The corresponding 
mean wind speed, U,^, for these two cases is 12 miles/hour 
and 30 miles/hour, respectively. Figure 19 confirms that 
the predominant influence of the change in Reynolds' number 
brought about by increasing the initial friction velocity, 
is confined to a region near the wall, with the effect being 
most significant at the top of the ellipse where the velocity 
gradients and shear stresses are large. The other profile 
is calculated at a position upstream near the forward 
separation region where the adverse pressure field causes 
the velocity gradients and stresses to be relatively smaller. 

The Eff e ct of Rou ghne ss and Ellip s e G eo metry on Separation 

The severe adverse pressure gradients which exist up- 
stream from the body cause separation of the boundary layer 
near the stagnation region. Since the numencaL solution, 
as discussed in Chapter lit, gives a reasonably accurate 

prediction of the point at which this separation occurs, a 
characteristic length of the forward separation bubble can 
be determined. This calculation has been carried out for a 
number of elliptical aspect ratios and surface roughnesses 
and is given in .Figure 20. Examining the effect of rough- 
ness on the length, of the separation region, it is 
apparent that a decrease m z q produces smaller separation 
regions tending to approach some asymptotic value for a 
smooth surface. This is expected since roughness reduces 
the velocity near the surface decreasing the momentum 
available to overcome the adverse pressure forces. 
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Figure 20. Elliptical geometry and surface roughness effects on the 
length of the separation regions. 




The effect of the aspect- ratio on t.lie Length of the 
separation is similar to that of the roughness, with the 
Longer ellipses producing smalLer separation bubbles. This 
fact can be explained by examining the pressure distribution 
over the cylinders as shown in Figure 14. page 65. As' 
pointed out in the discussion of that figure, the adverse 
pressure gradient upstream from the body is more severe for 
the smaller aspect ratio ellipse, leading to the conclusion 
that separation occurs further upstream for the shorter 
ellipses. f 

G eometrical Effects o n the Turbulent Viscosity |i 

It is not only important to determine the effects of 
the elliptical geometry on the mean velocity field, but: it 
is also desirable to examine the influence of the obstruc- 
tion on the turbulent characteristics of the flow; • 

This has been done in Figure 21 showing the develop- 
ment of the eddy viscosity profiles over the 10/1 ellipse. 
Starting far upstream the eddy viscosity profile, given by 
curve A, is calculated by Equation 15 from the initial 
logarithmic velocity distribution* Proceeding downstream 
under the influence of a mild adverse pressure gradient}, 
the next profile, given in Figure 21 by curve B, is at a 
position upstream of the forward stagnation point of about 
twice the ellipse height. If can be seen in comparison with 
curve A, that, the viscosity has decreased slightly in the 
lower region of the boundary Layer becoming larger near the 
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lopment of the eddy viscosity profiles over 





outer, edge. Curve C is calculated over the stagnation point 
where : :\a imp:. sed mean pressure gradient first becomes favor- 
able. This profile continues the- trend discussed above with 
the viscosity reaching a maximum near the outer edge. The 
profile given at point D, is calculated at a position down- 
stream of the maximum favorable pressure gradient and 
exhibits a drastic change from the pattern established by 
the first three curves, with the eddy values increasing to 
local maximum values near the wall and decreasing rapidly to 
minimum values in the outej boundary layer. And, finally, 
the last profile, calculate, at the top of the ellipse 
(point E) , well downstream from D, shows the effect of the 
very small pressure gradient existing over a large portion 
of the cylinder. The response of the eddy viscosity to the 
influence of this negligible pressure field is similar to 
that of the velocity, tending to return to the form of the 
undisturbed flow upstream. This fact is not surprising, 
however, because the eddy viscosity, calculated by the 
Prandtl mixing length model, is directly related to the 
vertical gradient of the velocity. 

Therefore, the results given in Figure 21 can be 
explained by examining the effect of the pressure forces on 
the velocity. Since adverse pressure decelerates the flow 
reducing the gradient in velocity near the wall while 
increasing it in the outer flow, the turbulent shear will 
show a similar effect as illustrated in comparing curves A, 

B and C of Figure 21. On the oth-;r hand, the favorable 
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pressure gradients over the front of the ellipse will have 
the opposite effect, as shown by curve D with the large 
increase in eddy viscosity occurring near the wall decreas- 
ing rapidly in the outer boundary layer. 

Calculation of Flow Over a Fence with Comparison to Data 

In addition to the results previously discussed for 
atmospheric flow over the elliptical cylinders, several 
calculations have been carried out for flow over a fence, or 
two-dimensional vertical flat plate, which represents the 
limiting case of the elliptical cylinder with zero aspect 
ratio. This particular flow situation has a number of 
practical applications as a wind-sheltering device and has, 
therefore, been the subject of a number of experimental 
studies (4) , (11) and (12) . 

Applying the concepts developed within this study to 

flow over the fence, approximations of some representative 
data (12) can be made. Results of this calculation are 
given in Figure 22 and are compared with measured velocities 
upstream and on top of the fence. All calculations are 
carried out over a smooth surface using the Van Driest— 
Clauser eddy viscosity model with pressure distribution and 
outer boundary conditions given by the potential flow values 
along the first nonseparating streamline. The particular 
potential solution used in these calculations is obtained 
from Reference (23) assuming a symmetric flow over the 
fence . 
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Examining the two profiles upstream of the fence it 
appears that the numerical solution overpredicts the 
velocity near the wall with smaller predicted values in the 
outer edge of the flow. This is probably due to the large 
pressure gradient physically existing in the vertical 

m 

direction which decelerates the flow near the surface to a 
much greater extent than the assumed mean pressure gradient 
along the nonseparating streamline. The calculated profile 
at the top of the fence, while exhibiting a similar form to 
the measured profile, predicts much larger velocities in the 
boundary layer. This effect is due to the outer boundary 
conditions obtained from potential theory which do not take 
into account the very large separation regions both upstream 
and downstream of the plate. Because of the large decelera- 
tions produced by these regions, a potential solution which 
includes these effects would be a better approximation of 

the actual flow conditions, and more accurately predict the 
magnitude of the velocity. 

Both of the above effects produced by the large 
separation regions which exist in a major part of the flow 
field surrounding the fence, can be taken into account to 
some extent by the methods introduced in Chapter II. Since 
the pressure gradient imposed on the boundary .layer by 
reversed flow regions is much larger for the ease of the 
fence than for the elliptical cylinder, the former flow 
situation provides a severe test for the methods proposed in 

this investigation. However, with improved approximations 
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of the pressure field and boundary con- ,ions, it is 
believed that the boundary layer approach can make reason- 
able predictions of the flow field over the fence. 
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CHAPTER V 


DISCUSSION OF RESULTS FOR AN IMPROVED BOUNDARY 
LAYER ANALYSIS OF FLOW OVER 
ELLIPTICAL OBSTRUCTIONS 

Analysis of Enlarged Flow Fields 

The consideration of flow fields in a larger region 
above the elliptical obstruction than that discussed in the 
preceding chapter requires approximations of the outer 
boundary condition on the mean wind and the pressure dis- 
tribution over the obstacle which are more characteristic 
of atmospheric flow in this enlarged region than conditions 
prescribed along the nonseparating streamline. As discussed 
in Chapter II, the outer boundary condition which is a 
realistic approximation of the flow far from the surface is 
given by an asymptotic matching of the velocity within the 
boundary layer to the logarithmic velocity associated with 
the undisturbed flow. In conjunction with this boundary 
condition it was also pointed out that the pressure gradients 
associated with the obstruction must decay to zero from the 
value given along the nonseparating streamline. Therefore, 
the initial approximation of the pressure field existing 
within the region of influence of the elliptical cylinder, 
is given by the following second order quadratic decay 
function : 
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where ^ is the pressure gradient determined along 
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the nonseparating streamline over the obstruction. 

Preliminary calculations of the flow field in an 
enlarged region over the elliptical cylinder, using the 
above pressure distribution and logarithmic outer boundary 
condition, have been carried out for flow over the 2/1 and 
4/1 aspect ratio ellipses with different roughnesses. The 
results of these calculations for the 2/1 ellipse are given 
in Figures 23 and 24 showing the development of the boundary 
layer from a position upstream of 10 times the ellipse 
height to the top of the cylinder. This initial location is 
chosen in order to reduce the computing time required for 
these preliminary calculations. The logarithmic boundary 
condition is applied at a constant height above the surface 
of five times the ellipse height, with the solutions non- 
dimensionalized by the logarithmic velocity defined for 
convenience at a position of three times the ellipse height. 
One observes from the curves that the flow is initially 
decelerated from the starting profile at x/b = -10.0 through 
the range of approximately x/b = -3.0 and then is accelerated 
to the top of the ellipse, x/b = 00. Similar results are 
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Figure 23. Development of the boundary layer over a 2/1 

ellipse with variable pressure gradient in the vertical 
direction and logarithmic outer boundary condition for 
z Q /b = .005. 
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Figure 24. Development of the boundary layer over a 2/1 

ellipse with variable pressure gradient in the vertical 
direction and logarithmic outer boundary condition 
for z Q /b = .020. 



shown in Figures 25 and 26 for the 4/1 ellipse. 

Comparing these results with the velocity profiles 
presented in the previous chapter, the most significant 
effect of the vertical variation of the pressure gradient is 
the pronounced maximum in the wind velocity produced locally 
near the surface for the calculations made over the top of 
the ellipse. This result is due to a combination of the 
large accelerations created by the favorable pressure 
gradients which exist in this region of the flow field, and 
the influence of the decaying pressure function imposed on 
these gradients which decreases the acceleration, allowing 
the flow to return to the logarithmic velocity at large 
heights above the surface. 

The acceleration of the boundary layer occurs from a 
position x/b = -3.0 to the top of the ellipse where the 
maximum wind velocities are attained at an elevation of 
about 0.3 times the ellipse height. In magnitude the wind 
is accelerated to about twice its initial value. This com- 
pares to calculations discussed in the previous chapter in 
which the maximum velocity at the top of the 2/1 ellipse was 
only 35 per cent greater than the initial wind speed. 

However, returning to the solutions given in Figures 
23 through 26, consider the flow far upstream. Since the 
pressure disturbances created by the obstruction are small, 
the effect of the vertical decay of the pressure gradient is 
negligible, with the form of the velocity profiles decreasing 
only slightly from the initial logarithmic distribution. 


89 



Oi o 










•/. 

f) 


r> . o 


2.0 


1.0 


.5 


.2 


.1 


.05 


.02 



! u U 

Figure 26. Development of the boundary layer over a 4/1 

ellipse with variable pressure gradient in the vertical 
direction and logarithmic outer boundary condition for 


z/b - 0.020. 
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Proceed! ruj downstream, the adverse pressure increases, 
reaching a maximum at approximately x/b = -3.0, where the 
deceleration of the flow near the surface causes the 
velocity to significantly deviate from the logarithmic 
distribution extending the influence of the body to very 
large heights. 

The results given in Figures 23 and 24, pages 87 and 
88, and Figures 25 and 26 are, in general, consistent with 
the wind profiles given in qualitative discussions of flow 
over surface obstructions (2) and (4) and with some experi- 
mental measurements of boundary layers over fences (11) and 
(12) , in that the velocity near the obstruction exceeds that 
of the undisturbed flow. Moreover, this result is predict- 
able from the standpoint of the conservation of mass, which 
requires that the acceleration of the flow at the top of 
the body be a maximum close to the surface where the stream- 
lines are the most severely compressed. 

Therefore, based on the above considerations and 
intuitive reasoning about the distribution of the pressure 
disturbances above the body, it is believed that the method 
used in obtaining the results given m this section is a 
justifiable approach in approximating the flow fields over 
surface obstructions. However, additional experimental work 
is needed to further verify this concept and to more 
accurately determine the variation of the pressure field in 
the vertical direction. 
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Boundary Layer Effects Produced by Separation Regions 

The improvements made in the boundary layer analysis 
by simulating the effects of the separation regions with 
stationary, inviscid vortices, discussed in Chapter II, have 
been preliminarily investigated. The calculations, given in 
Figures 27 and 28, have been carried out for flow over the 
2/1 ellipse for several surface roughnesses with these 
results compared to the solutions given in the preceding 
chapter for flow over the same body without attached 
vortices . 

The effects of the simulated separation region on the 
development of the boundary over the ellipse is illustrated 
in Figure 24, page 88, by comparing the velocity profiles 
calculated at different x-stations along the surface. The 
influence of the vortices on the velocity is to increase the 
thickness of the shear layer by decelerating the flow in the 
boundary layer both upstream and downstream of the forward 
separation region. For the calculation carried out over the 
vortex, this same effect occurs in the upper portion of the 
profile, but near the vortex, the additional momentum 
imparted to the flow by the slip velocity, u s , along the 
dividing streamline increases the velocities m the lower 
boundary layer. It is useful in examining the effect of the 
separation regions on the f Low to consider the potential 
velocity field over the elliptical cylinder given in Figure 
7, page 17. Comparing the distribution along the non- 
separating stream lines ( |f — j- r - 0.60) for the cases given with 
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Figure 28. Effect of vortices on velocity profile at top of 
2/1 ellipse with different roughnesses. 
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and without the effects of the potential vortices, it can be 
observed that, both upstream of the forward separation region 
and at the top of the ellipse, the vortex decreases the 
potential velocity. This effect can be interpreted through 

•the Bernoulli equation as an increase in the pressure gra- 
dient in these two regions of the flow. Since the pressure 
fields at these two positions are adverse, this fact leads 
to the conclusion that the velocity in the boundary layer 
will be less for calculations which include the influence of 
the separation regions. However, the velocity is affected 
to a greater extent upstream where the difference in the two 
potential solutions is greater than at the top of the 
ellipse . 

This effect is further illustrated in Figure 28 for 
calculations made at the top of the obstruction, for differ- 
ent surface roughnesses, where the maximum decrease in 
velocity brought about by* including the vortices in the 
potential solution, is only about two per cent. An addi- 
tional observation from this figure is that the separation 
regions appear to have no effect on the difference in 
velocity brought about by changes in roughness, that is, the 
increase in velocity produced by the change roughness from 
z Q = .010 to .005 meters is very nearly the same as the 
solution obtained with the addition of the vortices. 

From the above considerations it can be concluded 
that this method of simulating the additional effects of the 
separation regions with the inviscid vortices attacks the 
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most difficult area of the flow field to predict by boundary 
layer analysis and suggests that the effects of these regions 
on the overall flow field are small. In addition, steps 
have been taken to obtain a more realistic physical inter- 
pretation of the solutions over the separation bubble by 
including the slip velocity as- a lower boundary condition on 
the equations. 

It is believed that this initial attempt at incorpo- 
rating the effect of the separated flow into the boundary 
layer approach will improve the approximation of the flow 
field near the separation regions, but preliminary investi- 
gation of this method indicates that further calculations 
and experimental work are needed to improve the assumptions 
made. 
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CHAPTER VI 


DISCUSSION OF RESULTS RELATIVE TO 
AERONAUTICAL APPLICATION 

Winds and turbulence are significant factors influ- 
encing the safety, quality, and efficiency of air transpor- 
tation. Loads due to extreme turbulence or unexpected gusts 
of wind can exceed the static strength of the aircraft or 
upset the handling qualities to the point where control is 
lost and accidents occur. The effects of moderate turbulence 
although generally not as catastrophic, can degrade the 
precision of flight and tax the fatigue strength of the 
structure . 

Donely (31) reports that approximately 80 per cent of 
rough air experience occurs at altitudes below 25,000 ft. 
and is consequently prevalent to short haul or feeder line 
aircraft and to long haul aircraft during climb and descent. 
In the terminal areas where precision flight is most critical 
the air motion is frequently more severe due to wind dis- 
turbances from buildings and topological surface features. 

The wind gradient or shear can either add or absorb 
energy from an ascending or descending aircraft. Since the 
aircraft is essentially a constant speed machine, a gradient 
wind can thus increase or decrease the rate of descent and 
the angle of descent as well as the margin for the stall 
safety. These problems tend to increase as the landing 
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speed decreases, and for V/STOL vehicles, the angular varia- 
tion with altitude can be a significant factor in compli- 
cating the pilot's task of controlling the aircraft. 

Automatic landing systems are also susceptable to the 
effects of turbulence. McManus (32) reports on early 
experience with the Trident automatic landing system: "When 

the Trident was first landed automatically in moderate 
turbulence, several instances of poor performance occurred. 
These were caused by significant single horizontal gusts." 
"These gusts were not predicted by the general turbulence 
models as they were special to the site." In one case, a 
moderate crosswind was blowing over a factory into the under- 
shoot area. Close inspection showed a long hanger with a 
double apex roof at approximately 45 degrees to the wind, 
directly in line with the undershoot area. This generates a 
stable wind pattern which a single stationary measuring 
station would not record as high turbulence, but an aircraft 
passing through the pattern would be subject to large air- 
speed variations. 

The problems of an aircraft flying through a head 
wind that varies with elevation has been investigated in a 
paper by Etkin (33) . The equations of motion were set up 
and solved for a flight at constant air speed. The effect 
of this wind gradient on the flight path during descent and 
climb is shown. It is seen that the deviation from the case 
with no wind gradient is quite considerable. 

In the case of descending from a given height into a 
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positive wind gradient (increasing wind with height) , an 
airplane overshoots the initial tangent by a distance s, as 
shown in Figure 29. For a linear wind variation, this dis- 

t 

tance is given as: 

1 3w „ .2 . 

3 = 2 n V A 1 sln ♦ 


where 

dw 

= head wind gradient 

= Airplane Air Speed 
t - time 
<j> = path angle 

The problem of the climbing aircraft differs from the 
glide in that the angle <f> is negative and that there is an 
additional force, the thrust. The path is a parabolic as 
before, the convexity in the same sense, so that the hori- 
zontal distance to reach a given height is concentrated as 
compared with the case of zero wind gradient; see Figure 29. 
The effect of a negative wind gradient (decreasing wind with 

height) on glide and climb is shown in Figure 29. 

■ % 

Another effect arises, when the wind gradient is not 
coming directly from ahead, but is blowing at a certain 
angle 6 to the flight direction. It is shown in Figure 30 
that a cross wind induces a lift addition on one side and a 
lift subtraction on the other side of the wing, thus pro- 
ducing a roll-moment. As w(z) or w sin 3 is a function of 
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elevation, the roll-moment will also vary with height 
according to w(z). For strong gradients of w(z), corre- 
spondingly strong gradients in the roll-moment must be 
expected. 

An approximate quantitative estimate of the influence 
of gusts occurring due to the presence of buildings is 
possible from the results of this study shown in Figures 23 
through 26, pages 87, 88, 90, and 91. These figures show 
relatively large wind acceleration or gusts near the top of 
semi-elliptical surface protrusions. Taking the lift as a 
typical example, the following analysis for a 2:1 ellipse 
with the velocity profile shown in Figure 23 provides a 
quantitative feel for these effects. 

Lift is given by the expression: 



The ratio of the lift in the disturbance of the building, 
Lg, to that in the natural wind, , ali other parameters 
being constant, is: 



A plot of the ratio Lg/L N versus elevation at the x/b = 0 
station for the 2:1 ellipse with z Q /b = 0.005 assuming 
V A /U^ z=3 j = 1.0 is shown in Figure 31. The plot indicates 
that an aircraft flying into the wind over a long building 
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or surface obstruction of semi-elliptical configuration 
would experience an increase in lift exceeding 20 per cent of 
th6 value for wind over uniform terrain if passinq below a 
distance equal to the obstruction height z/b =1.0, (i.e., 

40 feet above the ground for a 20 foot obstruction) . 

Figure 32 shows the variation in lift ratio L_/L., 
with longitudinal position from the top of the ellipse to a 
dimensionless distance x/b = -10 upstream of the approaching 
wind where the wind is logarithmic. Parametric variation 
with z/b is given in the curve. 

An aircraft passing over the obstruction at a fixed 
height of z/b = 1.0 would experience an initial increase in 
lift of 20 per cent which would diminish and change to a 10 
per cent decrease at x/b = -3.0. Thus, a simple 15 foot 
high quansit hut located on the order of 600 feet from the 
end of a landing strip can generate a 30 per cent variation 
in lift in a distance less than 150 feet. The pilot's 
correction time at 50 mph landing speed is on the order of 
two seconds illustrating the hazards of unexpected wind pro- 
files due to surface projection. 

It is apparent from the previous discussion that a 
model for predicting wind velocity profiles about surface 
obstructions is important to the design of airplane control 
systems and to the design of airports themselves. The wind 
model presented in this report is a step in that direction 
since it provides transportation designers with insight to: 

1. Determining the safest distance for structures 
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from landing strips, at least, from the windward side. This 
is of particular significance to V/STOL ports in urban areas. 

2. Analyzing existing airports to warn of possible 
unexpected wind patterns and to define optimum locations of 
measuring stations to detect these wind patterns . 

3 . Improving mathematical gust models employed in 
automatic landing systems. Present gust model generally 
assumed a power law variation of the mean wind, for example 
Skelton (34) uses: 

U (z) = U( Z;L ) (z/z 1 ) P(z o ) 

A mathematical description for the flow patterns induced by 
surface obstruction is required to make these gust models 
more reliable. 

The use of the present boundary layer analysis of 
wind profiles is of course limited to a two-dimensional 
elliptical geometry and subject to the assumptions regarding 
eddy viscosity and pressure gradient decay. Moreover, the 
numerical nature of the solutions make them difficult to 
incorporate into mathematical gust models and into computa- 
tion of airplane characteristics. 

In view of the promising applications the. flow 
analysis over obstructions has, however, it is believed that 
additional effort to expand and improve the model is 
definitely desirable. The model should be: 

1. Extended to other geometries, 

2. Investigated using other eddy viscosity 
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correlations, for example, conservation of turbulence 
kinetic energy, 

3. Employed with other hypothesis regarding the 
pressure gradient surrounding the surface obstruction, and 

4. Compared to wind profile data from a simple 
experiment over an elliptical configuration. 
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CHAPTER VII 


CONCLUSIONS AND SUMMARY 

The boundary layer approach presented in this investi- 
gation provides a reasonable first order approximation of 
the atmospheric flow field which exists around buildings and 
other surface obstructions. This method has been applied to 
flow over elliptical cylinders, and, despite the lack of 
experimental data in this area, certain conclusions can be 
drawn about the relationship of the parameters involved in 
this problem. 

1. An increase in the elliptical aspect ratio 
decreases the wind speed within the boundary 
layer at the top of the ellipse and returns it to 
the logarithmic distribution characteristic of 
undisturbed flow . 

2. Increases in surface roughness affect the flow by 
decreasing the velocity in the boundary layer, 
with the most pronounced effect occurring near 
the surface of the smaller aspect ratio ellipse. 

3. The initial friction velocity which enters the 
boundary layer equations through the Reynolds 1 
number has a negligible effect on the overall 
flow for the range of Reynolds' numbers considered 
in this study. 
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4. 


A decrease in the elliptical aspect ratio and an 
increase in the surface roughness cause larger 
separation regions. 

5. The decay of the streamwise pressure gradient in 
the vertical direction produces localized maxi- 
mums in wind speed at the top of a surface 
obstruction, which are expected in physically 
real flow situations. 

6. The preliminary investigation of the simulation 
of separation regions with inviscid vortices 
indicates that this method appears to have only a 
small effect on the overall flow but does improve 
the approximation of the flow field near the 
stagnation regions. 

In summary, the method developed in this study 
utilizes the simplifying assumptions of boundary layer 
theory in an effort to analyze the complex flow field which 
exists in atmospheric motions over surface obstructions. It 
has been shown that reasonable approximations of the physi- 
cal conditions which exist in the atmospheric flow near the 
ground can be made without sacrificing the simplicity and 
savings in time and money realized with a boundary layer 
solution as opposed to a solution of the complete equations 
of motion. Finally, it is believed that, with additional 
experimental and analytical work to verify some of the 
proposed assumptions, flow fields over surface obstructions 
can be predicted with sufficient accuracy to make the 
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boundary layer approach a valuable tool in the analysis of 
atmospheric shear flows . 
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APPENDIXES 



APPENDIX A 


UNIFORM POTENTIAL FLOW PAST AN 
ELLIPTICAL CYLINDER 

The potential flow field around a two-dimensional 
elliptical cylinder can be obtained from the theory of 
complex variables. Consider a uniform, irrotational flow of 
an incompressible, inviscid fluid in the complex Z plane 
parallel to the X-axis with a velocity U. Defining the 
complex potential, w, as 

w (Z ) = <j> (X , Y ) + i\p (X , Y) 

where Z = X + iY, the velocity field can be obtained by 
differentiating the potential with respect to a coordinate 
direction. 


= JL& + ■ 3^ dw 3Z dw 

3X 3X 1 9X dZ 3X dZ 


Defining the complex velocity as the negative of this 
differential. 


v 


dw _ _ ^ 
dZ 3X 1 3X 


i v 


where u and v are the velocities in the X and Y directions, 
respectively. For the case of uniform flow parallel to the 
X-axis , 
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! 


v = U 


dw 

dZ 


or 


w(Z) = -UZ 


Consider a circular cylinder of radius ! R placed in' an 
arbitrary flow field of complex potential f(Z). From the 
circle theorem of Milne-Thomson (23) , the complex potential 
of the flow around the cylinder becomes 

_f R 2 ' 

w(Z) = f (Z ) + f If— 

where f is the complex conjugate of f. For a circular 
cylinder in a uniform flow field, the complex potential is 


w (Z ) 



(36) 


To transfer the region of flow in the Z plane out- 
side the circular c>_inder into a region outside an ellip- 
tical cylinder in the z plane, as shown in Figure Al, the 
Joukowski transformation (23) should be applied to the above 
equation. This transformation has the following form, 

2 - 2 + w 


where 
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U.O) 


z Plane 


Figure A1 . Geometry for Joukowski transformation of flow 
over a circular cylinder to flow over an elliptical 
cylinder . 
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2 _ 2 ,2 

c — a — b 


R — f b) 


Solving for the inverse transformation, 
Z = ^ ( z ± / z 2 ~ ^7) 


(37) 


and applying to Equation 36, the flow field over the cir- 
cular cylinder of radius R = (a + b)/2 is transformed into 
flow over an elliptical cylinder of major axis 2a aligned 
parallel to the x-axis and minor axis 2b aligned with the 
y-axis, as illustrated in Figure Al. 

Using the positive root of Equation 37, the complex 
potential in the z plane becomes 


w(z) = -i- uj (z + / z 2 _ c 2) + + 

2 l z + ^ - c* l 

w(z> = |»ca + b)(& * * 2 V*V C " i ,38> 


Defining an orthogonal, elliptical coordinate system, 
(£,ri), the z plane can be transformed into a more useful 
form. This coordinate transformation from the (x,y) vari- 
ables in the z plane to (£,n) in the elliptic plane is 
defined by 


z = x + iy = c cosh (£ + in) 


or 
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z = c(cosh E, cos n + sinh E, sin nl 


(39) 


where the lines of constant £ are cofocal ellipses with 
major axis parallel to the x-axis, and lines of constant n 
are orthogonal hyperbolas. In Figure A1 , page 106, the line 
E, = E, q defines an ellipse with major axis of length 2*a and 
minor axis of length 2.b. From Equation 39 

x = c cosh E, cos n (40) 

y = c sinh ^ sin n (41) 

therefore 


a = c cosh £ o = \4 2 - b 2 cosh E, 0 

(42) 

b = c sinh £ 0 =~\Za 2 - b 2 sinh C 0 

(43) 


From the definition of the hyperbolic functions. Equations 
42 and 43 can be combined as 

i 


a + b = c ( cosh 

^0 

+ sinh 

£; 0 ) = ce^o 

(44) 

a - b = c(cosh 

^0 

- sinh 

K 0 ) = ce - S° 

(45) 


Introducing Equations 39, 44, 45 and the hyperbolic identi- 
ties into Equation 38, the complex potential for flow over 
an elliptical cylinder becomes 


w(z) = jU(a + b) 


c cosh ( S+in ) 


+ ^c 2 cosh 2 U+in) - c 2 


ce 


c cosh(g+in) - "Vc 2 cosh 2 (£+in) - c i 


ce 


-s< 
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U (a + b) 


exp(C+in-C ) + e ^ ) 

o o 

w(z) = U(a + b) cosh'£ + in - £ Q ) 

w(z) = 4> 4- iip = u(a + b) cosh[(? - £ ) + in] (46) 

Expanding and equating the real and imaginary parts, the 

stream function will be defined as 


i|> = U (a + b) sinh(£ - £ Q ) sin n 

ip = U'a + b) [sinh £ cosh £ - cosh £ sinh £ ] sin n 

0 0 (47) 


Using Equations 42 and 43 and expanding cosh £ q as an 
exponential and hyperbolic sine term, the following relation 
is obtained for the stream function: 

\ p = -Ub 1 |4| e ^ sin n + ul4^ - b^ sinh £ sin n (48) 

Nondimensionalizing this equation by U and b and defining an 
ellipse aspect ratio k = a/b, the stream function becomes 


Ub 



k + 1 
k - 1 



sin n + 



l)sinh £ sin n 


(49) 


This relationship defines the streamlines for uniform 
potential flow over an elliptical cylinder. 

The velocity potential function, <J> , is similarly 
derived and is given by the following equation: 
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(50) 


4> _ k+1 -K 


Ub 


k - 1 


cos n + k - 1 cosh £ cos n 


Since the velocity in any direction, s, is given by - 
the potential velocity can be obtained by differentiating 
Equation 50 with respect to the direction of the streamline. 
However, it is more convenient to calculate the velocity 
along the streamline in terms of the orthogonal elliptic 
coordinates K,n). 

Consider an incremental distance, ds, along the 
streamline as shown in Figure A2 in terms of the Cartesian 
coordinates (x,y) . 


(ds) 2 = (dx) 2 + (dy) 2 


(51) 


Since under the elliptic transformation, 

x = x (£ ,n) 
y = yK,n) 

Equation 51 becomes 


(ds) 2 = d£ + dn I + 
ot, an ) 


fl£ d? + dn I 


H 


3n 


or 


(ds) 2 = E(d£;) 2 + F d£ dn + G(dn) 2 


(52) 


where E, F and G are scaling factors defined by 


f -wl Z I 

E Mf) + 1 


Hi 


( 53 ) 
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X 


x + dx 



Figure A2. The relationship of differential distances in 
Cartesian and elliptic coordinates. 
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p _ 9x 9x 9y 9y 

35 3n H 3n 


( 54 ) 


G = 



(55) 


For an orthogonal system, F = 0; therefore, 
(ds) 2 = E(d£) 2 + G (dn) 2 


Along a curve, n = n , dn •= 0. 

ds - /E d£ (n — n ) 

o 

Similarly, along E, = £ , d£ = 0 
ds = /g dn (e = £ 0 > 

For conformal transformations, the scaling cit any point must 
be the same in all directions, consequently, 

/E = /G 


and 


9_ 

9s 




o 


_1 

/E 


9 

9n 


1 1=£>. 



n=n 


o 


1 9 


/E 




n=n. 


(56) 


(57) 


From Equations 56 and 57 the velocity components in 
the positive E, and n directions as shown in Figure A3 are 
given by 
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( 58 ) 


„ - 3<Hg,n) 

5 " ' a ac 

U i. (59) 

n /E 

From Equations 50 and 53, these equations become 


u 


C _ 


u 




2 2 2 
£ cos n + cosh E, sin n ^ 


(VP4 


e ^ cos t\ 


2“ ^ 

k - 1 sinh E, cos n [ 


(60) 


u 


n _ 


U -i/ . , 2_ 2 ^ . 2 r . 2 , 

ysinh E cos n + cosh £ sin n v 




e ^ sin n 


+ V k^ - 1 cosh £ sin n] 


( 61 ) 


Equations 49, 59 and 61 determine the potential solution for 
flow over an elliptical cylinder. 
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APPENDIX B 


UNIFORM POTENTIAL FLOW PAST AN ELLIPTICAL 
CYLINDER WITH FIXED VORTICES 

Consider a circular vortex of radius A rotating in a 
stationary, inviscid, irrotational, incompressible fluid. 
Locating circles concentric with the vortex of radius r* < A 
and r" > A, Stokes circulation theorem (23) can be applied 
to obtain the velocity of the fluid on these circles. 


1 

*2 

ojA 


q 2 

r h 

r > A 

1 

q = 2 

rw 

r ' < A 


where u> is the vorticity of the fluid inside of the vortex. 
Therefore, the existence of a vortex in an inviscid fluid 
induces a velocity field in the region around it with the 
magnitude of the induced velocity inversely proportional to 
the radial distance from the center of the vortex and 
directed perpendicularly to this radius. The fluid within 
the vortex rotates as a rigid body with a velocity propor- 
tional to the radial distance. It is important to note that 
the vortex cannot induce motion at its center, implying that 
a circular vortex in an undisturbed fluid will remain fixed 
at its initial position. 

Considering a circular vortex centered at the origin 
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in the complex Z plane, the velocity field induced to a 
stationary fluid is given by 



Centering the vortex at an arbitrary point z Q and defining a 

1 2 

circulation strength k = jA u, being positive in the counter- 
clockwise direction, the complex potential becomes 

w (Z ) = iic Jln(Z - Z ) (62) 

o 

The rotational core of the vortex can be reduced to a 
singularity at the center by considering the concept of a 
vortex filament. Taking the limit as the radius of the 
vortex decreases to zero and the vorticity approaches 
infinity, such that the product A • uj remains constant, the 
circular vortex degenerates to a point with the external 
flow field still given by Equation 62. 

Consider a pair of vortex filaments of equal strengths 
but opposite rotations as shown in Figure A4 , located at 



Figure A4 . Stationary inviscid vortices behind a circular 
cylinder. 
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conjugate points Z^ and Z^ outside a circular cylinder 
|z| = A. The complex potential of the vortex combination in 
the absence of the cylinder is given by Equation 61 as the 
sum of the potentials due to each vortex rotating alone in 
an infinite, undisturbed fluid. If the motion of the fluid 
is due solely to the vortices, the circle theorem gives the 
potential for flow around the cylinder as 


w(Z) 


-iictn 



+ 


iicJln 




(63) 


Expanding the logarithms of Equation 63, the potential 
becomes 


w (Z ) 


i<Jln(Z - Z.) + i<£n(Z - Z,) + iK&n hr— 


- i<ln \z~ ~ Z 1 



or 


w(Z) = -i<An(Z - Z 1 ) + iicJln(Z - Z^) 


+ 


i<tn 


A 2 ] 

z - 9- 

- iiciln 

f^i) 

z l] 


l z il 


iicJln I Z 



(64) 


The terms of the above equation irepresent, respectively, the 
potential of the vortex filaments at A and B, the images of 
these vortices at the inverse points A' and B' inside the 
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cylinder, plus a constant term. 

Since a vortex cannot induce velocity at its center, 
the motion of a vortex system is a function of the combined 
velocity fields of all sources, sinks, and other vortices in 
a fluid region. From Equation 63, the complex velocity of 
the vortex at A can be calculated by considering all contri- 
butions to the flow except from the vortex. Therefore, 
eliminating the first term of Equation 64 the complex 
potential of the vortex at A is given by 



(65) 


Adding a uniform velocity field, U, directed in the 
positive X-direction, to this flow field around the cylinder, 
the complex potential at A becomes 


W A (Z) 



( 


+ i<An 



■ 




( 66 ) 


The complex velocity at A is then obtained from the above 
relation by differentiating with respect to Z and setting 
Z = z 1% 

To fix this vortex at point A, the velocity, v A , must 
vanish. Applying these conditions to Equation 66 and 
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expand i n<j 


U 


l l - 6 


= IK* 


(Z x - A 2 ) (Z^ - A 2 ) + A 2 (Z 1 - Z x ) 
(Z^ - A 2 ) - Z x ) (Z 2 - A) 


dividing both sides of this equation by the conjugate of 
Equation 65 and simplifying, the condition for a fixed 
vortex at A becomes 


(Z 1 Z 1 - A 2 ) + Z 1 Z 1 (Z 1 - Z x ) 2 = 0 


(67) 


ie , 


Expressing Equation 67 in polar form, Z^ = re it becomes 
(r 2 - A 2 ) 2 + r 2 (re 10 - re -10 ) 2 = 0 


2 , 2 


4 • 2 , . 


2 4 . 2, 


(r - A ) = -4r sinh ( i 0 ) = ~4i r sin 6 


2 2 2 
r - A = 2r* sin 6 


r - — = 2r sin 0 
r 


( 68 ) 


Referring to Figure A4, page 117, Equation 68 yields the 
following condition for a stationary vortex at A: 


AA* = AB 


(69) 


In addition. Equation 68 relates the strength of the vortex 
at A to its position by 


K A = 


U (r 2 - A 2 ) 2 (r 2 + A 2 ) 


(70) 
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Since the flow is symmetric to the X-axis, similar condi- 
tions apply to the vortex filament at B. 

The complex potential for flow of the uniform stream 
U over a circular cylinder with fixed vortices behind is 
obtained by adding the term -i<&n(Z - Z^) for the vortex at 
A to Equation 66, resulting in 


w(Z) 


-u|z + 



+ iicfcn (Z 


Z^) + iicin 



- i<£.n(Z - Z^) 


itcAn 



( 71 ) 


To include vortices in the flow field upstream of the 
circular cylinder at positions Z 2 and z" 2 , similar logarithmic 
terms must be added to the above equation yielding a complex 
potential of the following form: 


w (Z ) 


- -u(z + |ij + i-c 


£n(Z - Z^) + Jin 


Z - 


- Hn (Z 


- Z L ) - Jlnjz - f~] + *-n(Z - Z 2 ) 




,2i 


+ Jtnf Z - — 

- £n(Z - Z 2 ) - £n 

Z - §- 


1 


, z 2 l 

i 


( 71 ) 


Introducing the complex coordinates, Z = X + iY, into the 
above equation and separating the real and imaginary parts 
results in the following relation for the nondimens ional 
stream function: 
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* - 0A = Y - - *!>* + < Y - *i>“> 

A t x 


+ Jin 



r x - X 1 

2 

4 . 

Y i 

Y 

2 


2 ,,2 

T 




1*1 + 


x i + Y iJ 



JLn[(X - X 1 ) 2 + (Y - Y x ) 2 ] - In 


X. 


X + 


1 — 7 

x i +Y V 


Y + 


2 2 
X 2 + Y i; 


+ Hn [ (X - 


X 2 ) 2 + (Y + Y 2 ) 2 ] 


- Jin 


X - 


2 2 
X 2 + Y 2) 


Y - 


~2 2 
X 2 + Y 2 


+ Jin [ (X - X 2 ) 2 + (Y - Y 2 ) 2 ] 


+ Jin 



f 

V 

X 2 

CN 

CN 

H 

> 


2 


V _ 





r“ 

> 

1 

+ Y 2 
X 2 J 

T 

l ^ 

Y 2 

2 J 



) ( 72 ) 


Differentiating the stream function with respect to the X- 
and Y-directions provides the components of the velocity 
vector. 


u 


3Y 


( 73 ) 


v 



( 74 ) 


where 


These equations together with Equation 72 define the invis- 
cid streamlines and velocities over the circular cylinder 
with stationary vortices upstream and downstream in a 
uniform flow field. 

Applying the Joukowski transformation converts this 
region of the flow outside the circular cylinder in the Z 
plane to a region outside an elliptical cylinder in the Z 
plane. 


z = 



where 


c 


2 


b 


2 


In terms of the coordinate direction. 


+ iy = X + iY + 


2 ,2 

a - b 

4 (X + lY) 


Nondimensionalizing by the length of the minor axis 2 • b 
and defining an ellipse aspect ratio, k = a/b. 


x 


X + 



- 1) * X 

~ir — t~ 

IT + Y 


(75) 


y = Y - 


(k 2 - 1) • Y 

3 2 — 

X + Y 


(76) 


Applying the transformation given by Equations 75 and 
76 to the expression for the stream function of Equation 72 
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and the equations for the potential velocity of Equations 73 
and 74 gives the solution for inviscid flow over an ellip- 
tical cylinder with fixed vortices. 
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APPENDIX C 


THE CONSERVATION OF EDDY VISCOSITY 

A majority of the methods used for incorporating the 
effect of turbulence on the boundary layer have basically 
relied on the mixing length theory as an empirical relation- 
ship between the turbulent viscosity and the mean flow 
parameters. A more recent concept, proposed by Nee and 
Kovasznay in (21), is based to a greater extent on the 
physical mechanisms involved with turbulent interaction by 
considering the eddy viscosity as a property of the flow 
which is conserved during turbulent motion within the 
boundary layer. 

In describing the distribution of the turbulent shear 
through the boundary layer, the momentum equation relates 
the Reynolds' stress to the processes associated with con- 
vection, diffusion, production and dissipation. Assuming 
that the eddy viscosity is conserved during the exchanges of 
momentum involved with these processes, the following trans- 
port equation has been proposed 

- 3n . - 3n 3 (- 

3x 3z 3z [ 

n (n - v) QU °° , du 
U 2 dx 3z 

oo 

where the bar indicates dimensional quantities. The 


+ a (n - v ) 
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variable* it is the sum of the molecular and eddy viscosities; 
a, b and c are empirically determined but universal con- 
stants; and T* is a length scale associated with the dissi- 
pation of turbulence which is equivalent to the normal 
distance from a solid boundary for turbulent shear flows. 

Nondimensionalizing this equation results in the 
following nonlinear, partial differential equation for the 
eddy viscosity: 


u 


i£ 

3x 


+ w 


3e_ 

3z 


!_ 

Re 


(e + 1)1^ ♦ 


b 

Re 


3z 

e (e + 1) 


de 

3z 


+ ae 


c 

Re 


U 


3u 


3z 


e(e + 1) dU e 


dx 


3u 
3z I 


( 78 ) 


The boundary conditions which apply to this equation 
will, of course, vary with the flow situation considered, 
but for the normal boundary layer problem with an inviscid, 
turbulence-free external flow, the following conditions 
apply at the boundaries: 

e (x,0) = 0 

lim e (x, z) = 0 

y+oo 

The solution of this equation also requires an initial con- 
dition on the eddy viscosity. One possibility is to apply 
some type of mixing length model to an initial velocity 
profile . 

Numerical solution. Since the eddy conservation 


140 


equation and the momentum equation are coupled through the 
convection and the dissipation terms, respectively, the two 
solutions must be obtained simultaneously through an 
iterative technique. Applying an implicit finite difference 
method similar to the solution of the boundary layer equa- 
tion, results in a system of algebraic equations which can 
be solved by the inversion of the tridiagonal coefficient 
matrix . 

The conservation equation is linearized by the 
velocity profiles calculated along the previous x-station. 
Applying the above numerical solution, the eddy profile is 
used to re-enter the momentum equation and solve for the 
velocity, iterating on this procedure until the solutions 
for the velocity and eddy viscosity converge. 

An attempt has been made to incorporate the concept 
of conservation of eddy viscosity into the turbulent boundary 
layer model, but problems in obtaining a satisfactory initial 
eddy profile leading to convergence of the solution and 
insufficient time to investigate the possibilities of this 
method, have prevented any meaningful results or evaluations 
of this theory. 
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